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FINAL  REPORT 


on 

MOLECULAR  INTERACTIONS  WITH  MANY-BODY 
PERTURBATION  THEORY 

to 

AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
from 

BATTELLE  LABORATORIES 
C0LUM8US,  OHIO 


I.  INTRODUCTION 

In  a  wide  variety  of  Air  Force  applications,  highly  detailed  information 

f 1-31* 

about  atoms,  molecules,  and  their  interactions  is  required.'  '  This  infor¬ 
mation  is  necessary  in  problems  ranging  from  chemical  laser  development,  to  the 
detection  and  identification  of  rocket  plumes,  to  metal  clustering  and  aerosol 
formations,  and  even  to  nuclear  weapons  effects. 

Probably  the  most  crucial  component  needed  to  understand  molecular  reactions 
is  the  potential  energy  surfaces  that  serve  to  describe  the  attractions  among 
the  atoms  and  molecules. ^  However,  such  information  is  not  easy  to  obtain. 

A  certain  amount  of  information  about  the  molecular  forces  near  equilibrium  in 
a  bound  molecule  is  available  from  spectroscopy.  Some  information  about  the 
potential  energy  surface  even  in  the  absence  of  binding  can  be  provided  from 
crossed  molecular-beam  experiments.  But,  in  general,  potential  energy  surfaces 
are  not  amendable  to  experimental  determination.  Instead,  other  types  of  ex¬ 
perimental  observations  such  as  kinetics  experiments,  coupled  with  very  simple 


theoretical  models  for  a  surface,  are  used  to  infer  pieces  of  information  about 

the  parameters  of  the  model  such  as  what  the  activation  barrier  might  be. 

The  most  direct  approach  to  obtaining  detailed  information  about  a  potential 

energy  surface  is  offered  by  predictive,  ab  initio  quantum  mechanical  calculations. 

However,  to  make  it  feasible  to  calculate  accurate  energy  surfaces  for  molecules, 

much  better  and  more  computationally  efficient  methods  must  still  be  developed. 

One  such  approach,  namely  many-body  perturbation  theory  (MBPT ) and 

its  infinite-order  extensions  termed  coupled-cluster  methods  (CCM) 

offer  a  number  of  attractive  features  that  the  more  traditional  configuration 

(21 1 

interaction  approaches  lack/  '  Under  AFOSR  support  at  Battelle's  Columbus 

Laboratories,  very  efficient  computer  codes  to  perform  MBPT/CCM  calculations 

were  written  and  employed  for  the  first  time  in  large-scale  ab  initio  calculations 

(11  21 1 

of  potential  energy  surfaces/  '  The  successes  of  this  effort  have  been 
substantial.  These  include  the  determination  of  a  complete  force-field  for  the 
H.^0  molecule,  including  all  force-constants  through  fourth-order,  that  is  suf¬ 
ficiently  accurate  that  once  improved  experiments  were  carried  out  after  our 

calculations,  many  of  the  previously  accepted  values  for  the  force  constants 

(22) 

were  revised  to  be  more  consistent  with  our  predictions/  '  Also,  a  study  of 
the  binding  energies  of  the  molecules  BgHg^BHg,  HgBN^-’-BH^+NHj,  and  H^BCCH-BHg+CO 
was  made  that  predict  these  binding  energies  to  within  1  kcal/mole  of  the  accepted 
experiments  for  diborane  and  borane  carbonyl,  and  make  a  prediction  in  the  case 
of  borazane  in  the  absence  of  an  experiment. Earlier  experiments  which 
gave  much  higher  values  for  the  binding  energies  of  diborane  and  borane  carbonyl 
than  we  computed  are  now  completely  discounted.  Similar  successes  with  studies 
of  the  isomerication  energy  and  activation  barrier  of  HNOHCN*^^  and  CH^NC/CH-jCN, 
the  photo-dissociation  of  formaldehyde,^5^  arid  various  studies  of  methanol, 
methoxy,  and  the  formyl  radical attest  to  the  r@l lability  of  our  MBPT/CCM 


methods. 
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Building  upon  this  work  supported  by  the  AFOSR,  at  Battelle,  we  carried 

out  extensive  studies  of  the  potential  energy  surface  for  the  two  inelastic 
3  3 

collisions,  0(  P)  +  F^O  and  0(  P)  +  COj,  under  contract  to  the  Air  Force  Rocket 

Propulsion  Laboratory,  for  the  purpose  of  obtaining  vibrational  excitation  cross- 

(271 

sections  that  are  needed  in  actual  detection  devices.  ' 

Despite  the  many  successes  we  have  had,  there  are  still  categories  of  prob¬ 
lems  that  cannot  yet  be  attacked  by  MBPT/CCM.  These  include  studies  of  most 

excited  states,  reactions  that  break  multiple  bonds,  and  applications  to  various 

(211 

kinds  of  open-shell  molecules.'  '  To  satisfy  these  additional  requirements 
it  is  necessary  to  simultaneously  develop  the  formal  theory,  write  additional 
computer  programs,  and  continue  to  make  landmark  applications  of  our  developing 
quantum  mechanical  technology.  Although  in  many  cases  the  formal  theory  is 
less  dramatic  than  the  applications,  the  continual  extension  of  the  theory  has 
a  greater  impact  on  our  ability  to  calculate  accurate  energy  surfaces  for  what¬ 
ever  categories  of  problems  might  emerge  from  the  needs  of  the  Air  Force. 

Consistent  with  this  objective,  much  of  our  work  has  been  devoted  to  formal 
theory.  This  includes,  the  derivation  of  the  coupled-cluster  single  and  double 
excitation  model  (CCSD)^a),  optimization  of  orbitals  within  the  coupled-cluster 
framework,  '^9)  an(j  developing  additional  mathematical  techniques  to  efficiently 
solve  the  nonlinear  coupled-cluster  equations. Additional  applications  to  a 
variety  of  problems  have  also  been  accomplished. 

In  the  following.  Section  II  discusses  the  research  objectives  of  this 
grant,  and  summarizes  some  of  the  notable  accomplishments  made  in  the  past  three 
years  under  our  AFOSR  grant  at  Battelle's  Columbus  Laboratories. 

In  Section  III  we  discuss  some  of  the  research  directions  that  work  in 


MBPT/CCM  should  take  in  the  future. 
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Section  IV  lists  the  publications  and  presentations  which  were  supported 
by  our  three-year  grant  with  AFOSR. 

II.  REVIEW  OF  RESEARCH  ACCOMPLISHMENTS 

Three  years  ago,  when  this  grant  was  initiated,  we  proposed  the  following 
overall  objectives: 

1)  Develop  new,  more  accurate  and  more  efficient  ab  initio  quantum 
mechanical  methods  based  upon  MBPT  and  CCM  for  determining  molecular 
properties,  and  particularly,  potential  energy  surfaces  for  molecular 
interactions. 

2)  Implement  these  methods  in  highly  efficient,  transportable  computer 
codes,  to  enable  computations  on  potential  energy  surfaces  to  be  made 
on  an  almost  routine  basis. 

3)  Apply  these  techniques  to  a  variety  of  problems  that  are  of  interest 
to  AFOSR,  and  that  serve  to  establish  the  range  of  accuracy  for  MBPT 
and  CCM  methods. 

In  line  with  these  objectives  a  number  of  firsts  have  been  accomplished  in 
this  program.  Highlights  are  reported  in  the  annual  reports  of  1979  and  1980, 
but  we  will  attempt  to  summarize  some  of  these  accomplishments  together  with 
those  of  the  past  year  in  this  final  report.  The  achievements  in  this  project 
have  been  reported  in  20  published  papers  and  30  presentations  including 
invited  lectures  at  most  of  the  principal  international  quantum  chemistry 
meetings.  These  are  listed  in  Section  IV. 

Some  of  the  principal  accomplishments  over  the  three-years  of  this  grant 
include  the  following: 

A.  The  first  general  purpose  implementation  of  the  infinite-order  coupled- 
cluster  double  excitation  model  (CCD)  was  made  and  applied  to  a  number 
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of  problems. >22,26)  unlike  Cl  methods,  which  for  practical  reasons 
are  usually  restricted  to  single  and  double  excitations,  CCD  properly 
includes  the  important  effect  of  quadruple  excitations  into  a  molecular 
calculation. 

B.  The  first  complete  derivation  and  implementation  of  the  coupled- 
cluster  single  and  double  excitation  model  (CCSD)  has  been  made  and 
initial  applications  recently  submitted  for  publ ication. The  CCSD 
wavefunction,  eTl+1'2  [<},  >,  includes  all  effects  of  single  excitations 
and  the  disconnected  triple  excitation  terms  such  as  T,T„,  and  the 

1  t 
4 

quartic,  T^  terms.  This  model  will  be  used  as  the  framework  for  a 
new  approach  to  chemical  bonding  using  localized  orbitals  where, 
unlike  SCF  orbitals,  T^  will  not  necessarily  be  small.  Since  CCSD  is 
equivalent  to  full  Cl  for  the  chemically  important  problem  of  separated, 
non-interacting  electron  pair  bonds,  we  expect  this  model  to  be  a  very 
interesting  study  in  our  future  work. 

C.  Both  the  infinite-order  CCD  and  CCSD  models  have  been  extensively 

appl  ied^1"27)  in  their  truncated  fourth-order  form,  termed  DQ-MBPT(4)  and 

SDQ-MBPT(4),  when  S,  D,  and  Q  refer  to  all  single,  double,  and  qua¬ 
druple  excitation  MBPT  diagrams  that  occur  through  fourth-order  in 
the  energy.  These  models  have  been  justified  by  showing  excellent 
convergence  to  CCD  and  CCSD  in  all  cases  where  SCF  orbitals  are  used 
and  no  quasi -degeneracies  are  encountered^21^. 

D.  In  other  work  involving  the  theory,  we  have  considered  orbital  op¬ 

timization  within  the  coupled  cluster  doubles  model,  developing  a 
quadratically  convergent  scheme.  In  future  work,  we  intend  to 

optimize  the  orbitals  for  the  CCSD  model. 
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E.  We  have  studied  the  isomerization  reactions,  HNC-HCN,  L i MOL i CN, 

( 2  2 )  ( 24 ) 

and  BNC-*-BCt!,  CH^NC+CH^CN.  txDeriments  exist  for  the  first  and  last 
isomerization,  while  we  predict  the  other  values.  In  tne  first  case, 
our  result  of  15  +  2  kcal/mole  was  in  disagreement  with  a  recent  ex¬ 
periment  that  obtained  10  kcal/mol^P^.Je  published  a  oaper  claiming 
the  experiment  was  in  error^.23^ecently,  Hehre  in  unpublished  work 
has  performed  an  ion  cyclotron  resonance  experiment  which  obtains 
14.3  ±  2  kcal/mole  vindicating  our  prediction.  In  the  case  of  the 
isomerization  of  methyl  isocyanide  our  prediction  of  22.7  kcal/mole 

is  in  excellent  agreement  with  an  experimental  value  of  23.7  obtained 

{ 32 ) 

by  Pri tchard1 s  group  in  Canada. 

F.  We  have  made  a  very  thorough  study  of  the  three  lowest  electronic 

surfaces  (X^A1,  a  ^A",  and  A*A")  for  HNO  H+NO  this  past  year^P^HNO  is 

a  molecule  common  to  many  flame  and  plume  species.  The  emission  of 
1 

HN0(A  A")-*-HN0(X  A1)  is  responsible  for  the  observed  red  emission  in 
the  reaction  of  H(  S)  with  N0(  IT).  Geometries  for  all  these  states 
are  obtained  providing  excellent  agreement  with  experiment  (within 

O  1 

0.02A  and  2°)  for  the  XxA'  and  A^A"  states,  while  predicting  the 

3 

structure  of  HNO  in  the  a  A"  state  in  the  absence  of  experiment. 
Excitation  energies  for  the  two  excited  states  and  recombination 

V  33  ] 

rates  for  the  H+N0-4IN0  reactions  are  also  obtained.' 

G.  The  first  application  of  CCD  to  a  potential  energy  surface  was  our 
determination  of  the  quartic  force  field  for  22 i n  this  work  we 
reported  all  force  constants,  made  extensive  comparisons  with  Cl  and 
with  different  MBPT  models.  Our  predictions  for  some  of  the  higher- 
order  force  constants  differed  sufficiently  that  Hoy  and  Bunker  re¬ 
interpreted  the  infra-red  data  using  a  more  flexible  treatment  of  the 
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bending  mode.  In  almost  every  case  their  revised  force  constants  were 
in  better  agreement  with  our  prediction. 

H.  Two  thorough  studies  of  the  formyl  radical,  HCO,  and  its  potential 
energy  surface  for  the  unimolecular  decomposition  HCOH+CO,  transition 
states,  heats  of  formation,  geometries,  and  activation  barriers  have 
been  reported. 

I.  This  year,  we  had  the  opportunity  to  participate  in  a  very  informative 
comparison  between  MBPT/CCM  and  full  Cl.  Full  Cl  refers  to  including 
all  possible  n-tuple  excitations  within  a  basis  set  for  a  molecule. 

Because  of  the  enormous  number  of  configurations  generated  by  the  Cl 
method,  such  a  solution  is  impossible  for  anything  but  the  simplest 
problems.  Saxe,  Schaeffer,  and  Handy  in  a  notable  achievement  have 

/  -JC\ 

obtained  the  full  Cl  within  a  double-zeta  basis  set  for  H^O. 

This  calculation  used  over  250,000  configurations,  and  required  about 
six  hours  on  a  CDC  7600.  Using  our  CCD  modei  which  includes  double-, 
quadruple-,  and  most  of  the  effects  of  higher-even  ordered  excitations 
and  adding  in  the  fourth-order  single  and  triple  excitation  contributions, 
we  obtained  a  result  within  0.2  kcal/mole  of  the  full  Cl.  Our  calcu¬ 
lations,  however ,  requi red  only  30  seconds  on  a  CDC  7600.  Although  this 
is  an  ideal  system  for  our  methods,  it  is  an  indication  of  the  high 

f  -v  -  \ 

degree  of  efficiency  attainable  by  many-body  methods.'1'3'  (See  review 
paper  on  MBPT/CCM  written  for  Annual  Reviews  of  Physical  Chemistry  for 
a  discussion  [36].) 

J.  The  structure  and  thermochemistry  for  the  highly  unusual  psuedo  tri¬ 
halogen  free  radical  HIF  has  been  obtained.  In  experiments  of  Yuan 

( 37 1 

Lee  and  coworkers  ,  this  radical  is  found  to  be  bound  by  30  kcal/mole 
relative  to  IF  +  H,  but  its  structure  is  completely  unknown,  without 
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even  any  evidence  whether  the  molecule  is  HIF,  HFI,  or  IHF.  Also,  it 
is  not  known  experimentally  whether  or  not  it  is  linear.  We  have  pre¬ 
dicted  its  structure  using  MBPT/CCM  and  effective  potentials  to  account 
for  the  chemically  inert  electrons  in  I.  1  The  molecule  is  found  to 
be  bent,  having  a  bond  angle  of  137°  with  I  in  the  center.  The  IF  and 
HI  bond  lengths  are  stretched  very  slightly  (0.1-0. 2ft)  from  the  cor¬ 
responding  bond  lengths  of  the  diatomic  molecules,  IF  and  HI.  We 
predict  a  binding  energy  of  25  kcal/mole  for  decomposition  to  IF  +  H, 
in  very  good  agreement  with  experiment. 

K.  The  reduced  linear  equation  method  developed  for  solving  coupled  cluster 

equations  has  been  generalized  by  investigating  a  least-squares  solution 

PC) 

to  the  psuedo-1 inear  CCM  problem.  ' 

L.  Comparisons  between  standard  configuration  iteraction  based  methods 
like  MCSCF,  truncated  Cl,  and  full  Cl  with  MBPT/CCM  methods  for  the 
insertion  of  Be  into  H^  have  been  initiated.  Even  though  this  system 
has  some  very  severe  degeneracy  problems,  our  initial  results  suggest 
that  a  single  reference  function  CCD  result  can  describe  this  insertion 
process  adequately.  The  system  is  also  serving  as  a  vehicle  to  illus¬ 
trate  che  orbital  optimized  CCD  method. 

M.  The  decomposition  of  formaldehyde,  H^CO,  to  radical  and  molecular  pro- 

(25) 

ducts  and  its  rearrangement  to  hydroxycarbene  has  been  studied.  ' 

This  problem  is  of  substantial  experimental  interest  because  of  formal¬ 
dehyde's  prominence  in  combustion/plume  processes.  The  activation 
barriers  and  heats  of  reactions  have  been  obtained.  In  the  latter 
case,  agreement  with  experiment  is  within  ±  2  kcal/mole.  The  activation 
barrier  predictions  support  the  Cl  results  of  Goddard  and  Schaeffer^  ^ 
that  would  suggest  a  tunneling  mechanism  for  H^CCWi^+CO. 
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N.  A  series  of  detailed  comparisons  of  various  MBPT  models  with  CCD  for 

the  C,  N,  and  0  atoms  and  the  f^O,  NH^,  and  CH^  molecules  have  been  made 

this  year.'  '  These  comparisons,  plus  a  number  of  others  we  have 

made,  suggest  that  the  infinite-order  CCD  results  differ  insignificantly 

(211 

from  the  fourth-order  model,  DQ-MBPT(4),  for  most  normal  cases.'  ' 

This  supports  the  predictions  of  the  less  expensive  fourth-order  model 
for  larger  molecules. 

O.  The  first  all-electron  ab  initio  coupled-cluster  and  MBPT  calculations 
of  benzene  were  made  last  year.  This  work  demonstrates  that  a  molecule 
of  this  size  has  at  least  a  20  percent  error  in  its  correlation  energy 
due  to  the  neglect  of  Cl  type  quadruple  excitations.^*^  This  emphasizes 
the  importance  of  using  methods  like  MBPT/CCM  that  properly  include 

such  higher  order  excitation  effects  if  reliable  quantum  mechanical 
calculations  are  to  be  possible  for  larger  molecules. 

P.  An  investigation  of  the  long-standing  discrepancies  in  the  experimental 

binding  energies  of  the  rocket  fuel  components,  BBH^-^Hg,  BH^+CO-H^BCO, 

(141 

and  BH3+NH3->-H3BNH3  has  been  made  with  MBPT.  The  binding  energies  of 
the  first  two  are  in  exceptional  agreement  (  ~1  kcal/mole)  with  the 
lower  set  of  experimental  values,  discounting  the  larger  binding  energies 
obtained  from  other  experiments.  The  value  of  borazane  of  30  kcal/mole 
is  a  prediction  in  the  absence  of  any  experiment. 


III.  FUTURE  RESEARCH  DIRECTIONS 

Over  the  duration  of  this  grant,  MBPT/CCM  relative  to  a  single  reference 
function  has  been  shown  to  be  a  quite  efficient,  accurate  method  for  ab  initio 
calculations  of  molecular  properties.  Numerous  studies  of  dissociation 
energies,  molecular  geometries,  and  force  constants  have  attributed  to  this 
fact,  as  highlighted  in  the  previous  section. 

MBPT/CCM  relative  to  a  single  reference  function  is  not  universally  ap¬ 
plicable,  however.  The  difficulty  does  not  lie  with  MBPT/CCM  but  rather  the 
limitation  to  a  single  reference  function.  In  our  work  we  have  chosen  to  use 
an  unrestricted  Hartree-Fock  (UHF)  function  to  permit  the  description  of  open- 
shells  and  to  correctly  separate  a  molecule  into  open-shell  fragments.  This 
seems  to  work  reliably  in  many  open-shell  cases  of  high  multiplicity  and  when 
breaking  a  single  bond,  but  the  correct  dissociation  of  multiply  bended  species 
is  more  difficult  to  describe. 

This  fact  is  illustrated  by  our  studies  of  the  potential  curve  dis¬ 
cussed  in  Appendix  A.  The  ground  state  of  is  ^g+»  but  a  restricted 
Hartree-Fock  reference  function  cannot  separate  correctly.  A  UHF  function  can, 
but  it  suffers  from  spin-contamination  since  the  spin  quantum  number  is  not 
conserved.  This  is  particularly  bad  for  a  singlet  state,  since  all  higher 
multiplicities  will  contaminate  the  lower  state.  Consequently,  as  disso¬ 
ciates,  instead  of  a  unit  multiplicity,  the  multiplicity  is  much  higher  ("-3.5) 
attributing  to  the  high  degree  of  contamination.  This  contamination  causes 
an  incorrect  behavior  in  the  UHF  description  of  the  region  where  M2  starts  to 
dissociate  that  subsequent  correlation  corrections  with  MBPT/CCM  cannot  cor¬ 
rect.  To  che  contrary,  there  are  few  observed  problems  in  studying  potential 
energy  surfaces  for  such  open-shell  species  as  HCO,  HNO,  CH^O,  and  0(^P)+H20, 
with  UHF  +  MBPT/CCM,  since  all  are  higher  multiplicities-  and  do  not  suffer 


from  the  same  degree  of  spin  contamination. 

There  are  potentially  three  solutions  to  the  problem  of  dissociation  of 
multiply  bonded  molecules  like  Ng-  The  first,  and  simplest,  would  be  to  find 
a  way  to  annihilate  the  spin  contamination  from  UHF  +  MBPT/CCM.  If  a  conve¬ 
nient  prescription  could  be  developed,  this  may  well  give  a  pragmatic  solution 
of  high  utility. 

The  other  two  possible  solutions  really  amount  to  the  same  thing  since 
each  introduces  additional  correlation  corrections  into  the  calculations,  but 
proceed  from  two  different  directions.  Either  one  can  retain  the  single  refer¬ 
ence  function  and  attempt  to  introduce  the  additional  higher-order  correction 
required  to  provide  correct  dissociations  such  as  coupled  cluster  single, 
double,  triple,  etc.  excitations (CCSDT. ..) ,  or  one  could  employ  multiple  re¬ 
ference  functions  which  should  not  require  as  high  an  order  treatment  of  the 
remainder  of  the  corrections.  The  last  approach  is  the  most  pervasive,  but 
also  requires  the  most  development  to  obtain  a  theory  that  is  convenient  com¬ 
putationally  and  generally  applicable.  Also,  the  so-called  quasi -degenerate 
forms  of  the  mul ti -reference  approach  suffer  from  serious  problems  with  intruder 
states,  perhaps  suggesting  that  a  fixed  linear  combination  of  reference  func¬ 
tions  might  be  a  preferable  starting  point. 

The  second  approach  using  a  single  reference  but  a  high  degree  of  corre¬ 
lation  has  been  applied  by  us  to  the  very  different  problem  of  the  insertion 
of  Be  into  H2  (see  Appendix  B).  This  problem  has  the  near  degeneracy  of  the 

2s  and  2p  orbitals  of  Be  as  well  as  the  degeneracy  encountered  when  the 

2  2- 

bond  is  broken,  where  the  lc  and  lcr^  configurations  become  equally  impor¬ 
tant.  This  system  is  small  enough  that  we  were  able  to  calculate  the  full  Cl 
(i.e.  all  possible  configurations)  solution  which  is  the  best  possible  solution 
for  the  problem „  The  CCSD  results,  even  relative  to  a  single  reference  function 


that  is  far  from  dominant  in  the  full  Cl,  is  essentially  indistinguishable  from  the 
full  Cl.  This  study  emphasizes  the  flexibility  and  stability  of  the  infinite- 
order  coupled  cluster  model  even  when  a  high  degree  of  degeneracy  is  present. 

For  more  complicated  problems,  however,  perhaps  mul tipi  e-reference  approaches  will 
be  the  only  solution. 

A  new  development  which  will  occupy  us  within  the  future  year  is  the  consid¬ 
eration  of  particular  classes  of  non-SCF  reference  functions  in  MBPT/CCM.  SCF 
reference  functions  provide  a  convenient  dichotomy  of  the  energy  for  a  molecule 
into  the  SCF  part  and  the  electron  correlation  effect,  which  is  recovered  by  MBPT/ 
CCM.  For  most  problems  SCF  theory  (at  least  in  the  unrestricted  form)  offers  a 
good  unperturbed  approximation.  Using  such  a  function  also  offers  a  number  of 
simplifications  since  all  nonvanishing  correlation  corrections  are  exclusively 
of  two-electron  type.  However,  SCF  theory  has  the  deficiency  that  the  canonical 
SCF  solutions  are  delocalized  over  the  entire  molecule.  For  larger  molecules, 
it  would  be  useful  if  the  orbitals  were  more  local  zed  since  many  of  the  two- 
electron  intergrals,  (aS|y6)»  would  essentially  vanish  for  charge  distributions 
ctB  and  v6  that  are  adequately  separated,  as  well  as  if  the  two  orbitals  in  the 
distribution  are  in  different  regions  of  the  molecule.  The  reduction  in  the 
number  of  nonzero  integrals  for  a  problem  can  have  a  drastic  effect  on  the  compu¬ 
ter  time  required  for  the  perturbation  calculations  for  large  molecules  as  well 
as  offering  a  more  conceptually  appealing  separation  into  units  related  to  electron- 
pair  bonds. 

Using  non-SCF  reference  functions  raises  the  question  how  high  an  order  in 
perturbation  theory  is  required  to  recover  the  same  quality  of  result  as  in  SCF- 
based  calculational  methods.  Also,  it  is  necessary  to  add  to  the  computer  codes 
all  the  additional  one-electron  perturbation  that  using  and  SCF  reference  function 
eliminates.  In  the  coming  year  we  hope  to  begin  to  answer  some  of  these  questions. 
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Another  area  that  we  intend  to  start  to  investigate  is  the  theory 


for  the  prediction  of  excited-state  potential -energy  surfaces  by  using  equation- 
of-motion  and  related  techniques  built  upon  a  MBPT/CCM  ground  state  reference 
function.  This  is  a  very  new  direction  for  this  project,  but  one  that  offers 
excellent  prospects  for  providing  excited  state  surfaces,  which  are  frequently 
important  in  classes  of  plume  and  combustion  problems.  As  a  bonus,  such  an 
approach  should  eventually  enable  us  to  obtain  electronic  excitation  spectra. 

The  fundamental  idea  is  that  for  the  wavefunction 


0 


V 


we  consider  a  second  operator,  r.,  such  that 


=  *1 


where  o-j  is  seme  excited  state.  Another  cluster  ocerator  is  a  possibility 
for,  such  as  .1  =  e^.  From  the  Schrodinger  equation, 


and 


H’<i  =  E  ’li 

o  o  o 


Hn  *  T 


so 


HC’i  =  £■  T.p 
0  1  0 


Left  multiplying  the  first  equation  by  n,  we  have 


[H,q]  =  aEo*0 


The  CC  wavefunction  yQ  may  be  obtained  in  the  usual  way,  while  a  set  of 
eouations  for  n  may  be  derived  from  the  equation-of -motion. 


In  addition  to  the  development  of  the  theory  along  the  lines  indicated 
above,  it  is  still  important  to  pursue  some  applications  to  problems  illus¬ 
trating  the  theory  at  different  levels  with  particular  regard  for  .-'■miparisons 
with  Cl,  full  Cl  when  possible,  and  some  MCSCF  results.  We  are  .  iterested  in  a 
variety  of  systems  like  interhalogen  flame  species.  In  particular,  the  unusual 
molecule  HCF  is  found  to  be  a  chemiluminescence  product  of  combustion  involving 
combinations  of  interhalogen  and  hydrocarbon  fuels,  as  may  be  important  in  rock¬ 
et  plumes.  We  intend  to  investigate  the  ground  state  energy  surface  for  this 
molecule.  Other  categories  of  molecules  to  be  studied  will  be  largely  determined 
by  examples  that  are  useful  to  illustrate  the  theory. 

We  intend  to  continue  investigating  these  possibilities  in  our  future  stud¬ 
ies  at  the  University  of  Florida,  Quantum  Theory  Project. 
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Abstract 

.1  lolecular  applications  of  coupled  cluster  and  manv-body  perturbation 
methods.  Rodney  J  Bartlett  and  George  D  Purvis  III  i  Bat  telle  Columbus 
Laboratories,  505  king  Avenue.  Columbus.  Ohio  43201.  L'.S.A.). 

Physica  Scripta  /Sweden)  71.  255-265.  1980. 

A  series  of  molecular  applications  of  manv-body  perturbation  theory 
I  M  RIM )  and  the  coupled-cluster  doubles  iCCDi  model  are  described. 
Even  rhough  these  methods  have  been  available  for  sometime,  only 
recently  have  large  scale,  MBPT  molecular  calculations  become  available. 
In  the  case  of  CCD.  the  results  presented  here  are  among  the  first  ob¬ 
tained  from  a  general  purpose  ab  initio  program.  The  intention  of  this 
paper  is  to  present  an  overview  of  the  current  state  of  the  manv-body 
approach  to  ground  state  properties  of  molecules.  The  properties  studied 
are  correlation  energies,  including  contributions  from  single,  double,  and 
quadruple  excitations  diagrams  in  fourth-and  higher-order;  dissociation 
energies;  potential  energy  surfaces;  and  molecular  polarizabilities  and 
hyperpolarizabilities.  Examples  are  taken  from  studies  of  a  variety  of 
molecules  including  HE.  H.O,  HCO.  C,fl(.  B;  H, .  CO.,  and  N. .  In 
many  cases,  it  is  found  that  quantitatively  accurate  dissociation  energies, 
geometries,  and  force  constants  can  be  obtained.  In  an  illustration  of  the 
•Y1 1  j  potential  energy  curve  of  N. .  it  is  shown  that  a  single  UHL  or  RUE 
reference  function  MBPT  CCD  approach  is  inadequate  at  some  inter- 
nuclear  separation. 

1.  Introduction 

Manv-body  (diagrammatic)  perturbation  theory  (MBPT)  basi¬ 
cally  dates  back  to  Brueckner’s  papers  [  1 .  2  J  of  1655  and 
Goidstone’s  proof  of  the  linked-diagram  theorem  in  1657  |3). 
Its  roots  include  the  work  of  Moeller  and  Plesset  in  1634  [4|. 
and  ultimately,  of  course.  Rayleigh-Schrodinger  perturbation 
theory  [5|.  The  applications  to  atoms  originated  in  Kelly's 
work  in  1663  [6| .  with  subsequent  applications  by  Das  and  co¬ 
workers  [  7  j  . 

Coupled  cluster  methods  (CCM),  developed  by  Coester  and 
Ktimmel  [8,  9]  in  1958-  1960.  and  in  a  form  useful  for  quan¬ 
tum  chemistry  by  Ci2ek.  Paldus.  and  co-workers  [10-131  can 
be  viewed  as  a  closed-form  set  of  equations  which  may  be  used 
to  sum  certain  categories  of  manv-body  diagrams  to  all  orders 
[14]  This  has  the  advantage  that  order  dependence  is  removed 
from  the  computation,  and  that  certain  invariance  properties  are 
present  that  would  not  normally  apply  to  a  finite-order  method. 
(Within  the  standard  MBPT  framework,  infinite-order  sum¬ 
mations  of  parts  of  diagrams  are  also  frequently  summed  to  all 
orders,  primarily  by  denominator  modifications,  but  not 
generally  entire  categories  of  diagrams.  See.  however,  (15)  and 
[16], I 

The  significance  of  MBPT/CCM  and  the  linked-diagram 
theorem  for  chemistry  has  several  facets.  Primary  among  these  is 
the  concept  borrowed  from  thermodynamics  of  "size- 
extensivity  '  [14]  The  term  indicates  the  proper  dependence 
>f  the  energy  or  density  matrix  on  the  size  of  a  homogeneous 


system,  and  is  a  necessary  result  of  the  exclusion  of  "unlinked" 
diagrams.  A  consequence  of  a  size-extensive  model  is  that  lor  a 
group  of  V  noninteracting  H;  molecules.  f.'(.VH:).  =  AVflH-) 
and  p(AH;)  =  .Vp(  H:).  Since  in  a  first  approximation  j  com¬ 
plicated  molecule  can  be  viewed  as  a  group  of  approximately 
nonimeracting  electron  pair  bonds  like  in  H: .  it  is  apparent  that 
the  size-extensive  property  should  be  maintained  as  the  theory 
is  developed  for  larger  and  larger  molecules.  A  truncated  Cl 
such  as  SD-CI  (all  single  and  double  excitations  from  a  reference 
determinant)  does  not  have  this  property  since  it  retains  un¬ 
linked  diagram  contributions.  For  SD-CI  i'(.VH;l  varies  as  the 
x/.V. 

Another  consequence  of  size-extensivity  for  chemistry  is  that 
for  a  reaction  consisting  of  closed-shell  species.  A’B-C-  D. 
the  heat  of  the  reaction  is  given  by  7sHrxn  =  AH,(C)  ~±H, 
(D)  --  \Hf{  A)  —  TsHfiB)  This  seems  like  an  almost  trivial 
result,  however,  if  the  heats  of  formation  of  the  species  were 
obtained  by  SD-CI.  this  simple  addition  is  not  entirely  justified. 
In  the  SD-CI  case,  one  would  prefer  to  obtain  the  heat  of  re¬ 
action  by  perUitming  "super-molecule"  calculations  of  the  two 
sides  of  the  reaction  with  A  and  B  and  C  and  D  infinitely  tar 
apart  to  partially  account  for  the  size-inextensivity .  If  A//f  is 
obtained  by  MBPT  CCM.  however,  a  table  of  theoretical  results 
for  individual  species  may  be  used  just  as  the  experimental 
values  are  employed. 

Although  superficially  similar,  the  correct  dependence  of  the 
energy  on  the  size  of  a  system  is  a  different  property  than 
correct  separation  of  a  molecule  into  its  fragments.  The  latter 
property  has  sometimes  been  called  "size -consistency  "  by  Pople 
and  co-workers  [I-.  18],  That  is  foi  any  molecule  AB.  a 
method  is  said  to  be  size-consistent  if  the  predicted  energies 
satisfv  f(AB)  =  E ( A)  +  E ( B).  where  A  and  B  max  be  open- 

R-» 

or  closed  shell  species.  It  is  apparent  that  if  the  linked-diagram 
expansion  is  not  truncated  or.  equivalently  if  a  fui!-(T  calcu¬ 
lation  is  made,  then  even  with  a  single  determinant  reference 
function,  size-consistency  is  guaranteed  The  ambiguity  between 
"size-consistency"  defined  as  correct  separation  and  "size- 
extensivity".  occurs  when  a  truncated  expansion  is  employed,  as 
is  necessary  in  any  practical  method 

The  truncation  of  the  expansion  toge  her  with  the  assump¬ 
tion  ot  a  single  determinant  reference  function,  forces  a  dis¬ 
tinction  between  closed-shell  molecules  separating  into  closed- 
shell  fragments  and  closed-shell  molecules  separating  into 
open-shell  fragments,  due  to  the  nature  of  the  restricted 
Uartree  I  t>ck  (RHFl  reference  function  (or  any  single  deter¬ 
minant  composed  of  doubly -occupied  orbitals  that  reflects  the 
point  group  symmetry  of  the  molecule  i. 

If  closed-shell  molecules  separate  into  closed-shell  f:  jcmettts. 
an  RI1F  Hincnon  is  a  proper  reference  function  for  jj!  utter- 
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nuclear  separations.  For  this  case,  any  approximation  to  ihe 
linked-diagram  expansion  for  the  energy  that  evaluates  an 
entire  diagram  (i.e..  not  just  the  "diagonal"  part  of  the  diagram, 
tor  example  because  of  possible  invariance  problems  [26]  )  is 
size-extensive  and  sue  consistent.  Hence.  EhrjH;  i  =  p/flH-  ). 

=  2£<Mg),  and  EJBjH*)  =  2£(BH3),  etc.  (Note,  the 
definition  ot  sue -consistency  or  correct  separation  only  pertains 
to  the  energy,  not  the  wavefuncuon  where  simultaneous  double 
excitations  on  the  various  fragments  are  necessary  even  in  the 
closed-shell  case.) 

However,  when  separation  of  the  closed-shell  molecule  into 
open-shell  fragments  is  required,  such  as  in  N2  or  H2,  the  RHF 
reference  tunction  for  the  molecule  does  not  change  smoothly 
into  a  single  determinant  RHF  functton  for  each  of  the  frag¬ 
ments.  but.  instead,  will  normally  go  to  an  ionic  form.  A*  and 
B~  Consequently,  there  is  no  single  determinant  RHF  reference 
function  to  employ  for  the  open-shell  fragment  in  a  consistent 
MBPT/CCM  calculation  to  investigate  whether  f(AB)  =  E( A)  f- 
F(B).  If  the  level  of  the  truncated  .MBPT  CCM  expansion  is  high 
enough,  even  though  the  RHF  reference  function  is  separating 
incorrectly,  the  correlation  corrections  are  sometimes  sufficient 
to  still  provide  a  good  potential  curve  at  large  imermolecular 
separations  [32] .  but  more  frequently,  it  is  necessary  to  require 
that  the  reference  function  should  also  separate  correctly. 

To  retain  the  simplicity  of  a  single  determinant  reference 
function  for  separation  into  open-shell  fragments,  one  must 
normally  resort  to  the  (spin  and  spatially)  unrestricted  Hartree- 
Fock  (UHF)  type  reference  function.  (Sometimes  different 
possible  UHF  functions  will  converge  to  different  separated 
atom  limits,  however,  so  care  should  be  exercised.)  Since  the 
UHF  function  will  normally  converge  to  an  RHF  function  for 
closed-shell  molecules  near  equilibrium,  the  UHF  function  will 
often  provide  a  reasonably  smooth  reference  determinant  as  a 
function  of  R  (see  Section  5  for  a  contrary  example  I.  but  once 
a  bifurcation  into  separate  RHF  and  UHF  functions  occurs,  a 
linked-diagram  expansion  can  be  evaluated  with  two  reference 
functions.  By  definition,  both  of  these  calculations  are  size- 
extensive  since  only  linked-diagrams  are  evaluated,  but  only 
one  of  these  two  calculations  would  permit  correct  separation, 
or  size-consistency.  Hence,  a  single  determinant  reference 
MBPT/CCM  calculation  normally  requires  a  UHF  function  to 
be  both  size -extensive  and  size-consistent,  and  heats  of  reaction 
with  open-shell  components  can  be  computed  accordingly. 

A  more  thorough  solution  to  the  ambiguity  between  correct 
separation  and  the  proper  dependence  on  the  size  of  a  system, 
requires  an  open-shell  MBPT/CCM  approach  [36-38] ,  since  a 
multi-determinant  reference  space  is  usually  required  to 
guarantee  correct  separation  in  a  truncated  expansion. 

The  property  of  size-extensivity  in  a  theoretical  model,  then, 
is  simply  a  consequence  of  a  more  proper  treatment  of  quad¬ 
ruple  and  higher  Cl  excitations  in  molecular  applications  [14, 
19.20].  These  excitations  are  responsible  for  the  cancellation 
of  unlinked  diagrams.  Thus  a  statement  that  size-extensivity  is 
important  is  simply  a  statement  that  quadruple  (predominantly) 
and  higher  Cl  excitations  are  important. 

These  higher-excitation  effects  are  handled  in  two  stages  in 
many-body  methods  [14.  19] .  The  first  stage  consists  of  incor¬ 
porating  higher-excitation  effects  to  eliminate  the  unliked  dia¬ 
grams  in  the  theory.  Thus,  any  approximation  to  the  linked- 
diagram  theorem  benefits  implicitv  from  this  feature.  The 
second  stage  comes  where  higher  order  Cl  excitations  like  quad¬ 
ruple  type  dV't’oh  are  further  decomposed  into  a  more  physi¬ 


cally  satisfying  set  of  components,  particularly  !,2  7";  <90' 
which  represent  two-simultaneous  double  excitations  [20]  By 
establishing  that  this  term  out  of  five  constitutes  the  pre¬ 
dominant  quadruple  excitation  contribution,  for  closed-shed 
molecules,  far  more  tractable  methods  for  including  most  -if  the 
effects  of  quadruple  excitations  are  possible  than  within  the 
standard  Cl  framework  [14] . 

In  Cl.  multi-reference  function  techniques  that  would  in¬ 
clude  all  single  and  double  excitations  out  of  the  reference 
determinants  would  introduce  excitations  which  are  four-fold 
relative  to  a  single  determinant,  and  as  such,  would  presumably 
contain  the  most  important  higher  excitation  effects.  Such  a 
method,  although  not  rigorously  size-extensive,  probably  would 
have  this  property  to  a  high-degree  of  accuracy . 

From  the  foregoing,  it  is  apparent  that  MBPT/CCM  has  mucn 
to  offer  in  molecular  problems.  However,  to  exploit  this  fact  in 
practice,  it  is  presently  necessary  to  use  conventional  finite  basis 
sets  of  Slater  type  orbitals  (STO)  or  contracted  Gaussian  type 
(CGTO).  This  is  a  consequence  of  the  muiticenter  nature  of 
molecular  charge  distributions.  The  few  attempts  to  use  one- 
center  expansion  techniques  for  molecules  containing  a  heavy' 
center,  which  permit  numerical  calculations  as  in  atoms,  met 
with  very  limited  success  [21.22].  Hence,  for  more  general 
molecular  environments,  basis  set  expansions  centered  at  various 
atoms  in  a  molecule  still  remain  necessary.  This  MBPT  CCM 
approach  has  been  pursued  by  Robo  [23] .  Bartlett  and  co¬ 
workers  [14.  15.  19.  20.  24.  25-31] ,  and  now  bv  several  groups 
[  1  *,  32-34] .  The  use  of  basis  sets  introduces  an  inherent  error 
in  the  calculations,  but  an  error  that  all  practical  quantum 
chemical  methods  share.  The  ultimate  answer  obtainable  in  a 
basis  set  is  the  ”fuU"  Cl.  Thus,  the  goal  -if  any  basis  set  method 
is  to  approach  as  closely  as  possible  to  this  result.  The  advan¬ 
tages  of  MBPT  CCM  for  including  important  effects  of  higher 
excitations  suggest  that  MBPT 'CCM  has  promise  of  converging 
more  quickly  toward  the  full  Cl  result  than  other  techniques. 

In  the  following,  we  attempt  to  prc  ide  an  overview  of  the 
current  state  of  MBPT  CCM  studies  of  molecular  systems.  Many 
of  the  results  reported  are  new  while  some  particularly 
illuminating  applications  made  by  our  group  have  jppearea  else¬ 
where  [14,  25.  28.  29]  After  a  brief  discussion  of  MBPT  CCM 
in  Section  2.  in  Section  3.  we  study  the  extent  of  the  corre¬ 
lation  energy  obtainable  in  a  standard  quality  basis  set  and  the 
contribution  from  the  different  orders  of  perturbation  theory 
Emphasis  is  placed  on  the  importance  of  quadruple  excitations. 
In  Section  4  the  prediction  of  dissociation  energies  with  MBPT 
is  discussed.  Section  5  describes  MBPT  CCM  applications  to 
potential  energy  surfaces,  with  some  comparisons  with  the 
SD-CI  model.  Section  6  describes  MBPT  applications  to  other 
properties  than  the  energy.  This  includes  results  from  the  first 
correlated  study  of  molecular  hyperpolarizabilities  [35] . 

2.  Summary  of  many-body  methods 

A  fairly  complete  discussion  of  MBPT  CCM  and  their  relation¬ 
ship  to  Cl  has  been  given  elsewhere  [14],  Consequently  ,  the 
present  section  will  be  limited  to  only  a  few  basic  equations  in 
order  to  define  some  terms. 

In  MBPT  the  energv  is  eiven  bv  the  linked-diaeram  expansion 

[3]- 

AE  =  E-E0  =  Et  = 

=  V  ot>oir[(£0  -  //or^Ti-Vr.  (11 
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/•Vi'.  I  An  nsy  mine  trued  (»<U)stune  diagrams  iAS(IDs)  through  fourth 
order  MBIT.  A  Hartree-I  oek  reference  state  is  assumed.  Orders  are  dis¬ 
tinguished  by  the  number  of  dashed  horizontal  interaction  lines.  Particle 
IP!  and  hole  ill)  states  are  represented  by  downward  and  upward  directed 
line  segments,  respectively  The  excitation  level  of  a  diagram  is  dis¬ 
tinguished  by  the  number  of  p-h  pairs  intersected  by  an  imaginary 
central  horizontal  line.  In  this  manner  the  diagrams  have  been  labeled  and 
counted  as  contributions  to  the  correlation  energy  arising  from  single  <S). 
double  i  D  l.  triple  1T1,  and  quadruple  (Q)  excitations  types. 


For  the  purpose  of  this  paper.  >b0  is  the  single  determinant  SCF 
result  for  a  nondegenerate  ground  state. 

The  Flamiltoman.  H0.  is  the  sum  of  one-electron  Fock 
operators  and  E0  is  the  sum  of  the  SCF  orbital  energies,  so  that 
t'ser  -  Ea  +  E, .  The  perturbation  is  V  -  H  —  Hn.  where  H  is 
the  usual  electrostatic  Hamiltonian,  and  the  subscript  L  indi¬ 
cates  the  limitation  to  linked  diagrams.  All  the  diagrams  (assum¬ 
ing  antisymmetrized  vertice  .)  that  need  to  be  considered  through 
fourth-order  are  shown  in  Fig.  1.  For  purposes  of  discussion, 
these  diagrams  are  characterized  by  the  number  of  particle  lines 
occurring  in  intermediate  vertices  into  single,  double,  triple,  and 
quadruple  excitation  types.  The  models  are  thus  defined  in  the 
form  SDQ-MBPT  [4]  which  means  all  single,  double,  and 
quadruple  excitation  diagrams  through  fourth  order  in  the 
energy  .  From  the  onset,  finite  basis  sets  are  assumed,  hence  the 
SCF  equations  are  understood  to  be  valid  in  a  matrix  sense. 

For  open-shell  atoms  and  molecules.  <!>„  is  chosen  to  be  the 
unrestricted  Hartree-Fock  (L'HFp-SCF  result,  where  relaxation 
of  spin  and  spatial  symmetry  is  permitted.  The  UI1F  solution 
will  normally  converge  to  the  restricted  (RHF)  solution,  which 
maintains  spin  and  spatial  symmetry,  for  gound  state,  closed- 
shell  molecules,  and  away  from  the  separated  atom  limit.  For 
many  examples  RHF  or  UHF  functions  provide  a  reasonable 
starting  point  for  a  correlated  method,  but  when  this  is  not 
possible,  such  as  when  more  than  one  determinant  would  be 
heavily  weighted,  it  is  necessary  to  use  multireference  function 
based  methods  [3b- 38)  .  In  our  work,  the  spin-multiplicity  for 
a  UHF  based  correlation  calculation  is  monitored  to  provide 
evidence  that  the  appropriate  spin  state  is  being  described. 
Although  it  is  better  to  treat  open-shell  problems  with  multi¬ 
reference  function  techniques,  the  greater  complications  in¬ 
volved  in  these  methods  makes  the  UHF  approach  an  attractive 
intermediate  level  technique  that  though  not  of  universal 
applicability,  is  of  wide  applicability. 

The  Rayieigh-Schrodinger  perturbation  theory  (RSPT)  ex¬ 
pression  for  is  [5 1  . 


2k£  =  V  <'!>„, K|(£0  -//„)' 1  PlV -Ah'  i|n  ;■!»„>  <2| 

P  is  the  projector  onto  the  space  orthogonal  to  'P,, .which,  in 
this  case,  represents  the  space  of  all  Cl  excitations  from  'l>(J 
Equation  (2)  is  equal  order-by-ordet  with  the  linked-dtagram 
expansion,  eq.  (I  ).  if  the  Cl  excitation  space  is  riot  truncated. 
However,  if  P  is  represented  b-z  only  double-excitations  from 
‘ho.  then  eq.  (2)  will  converge  to  the  D-(T  solution  which  re¬ 
tains  unlinked  diagrams.  This  should  be  contrasted  with  eq.  1 1  I. 
where  the  cancellations  of  all  unlmked-diagrams  has  already 
been  achieved  by  permitting  mixing  between  double  and  higher 
excitations.  This  cancellation  is  responsible  tor  the  sue  ex- 
tensivitv  and  the  particular  utility  ot  the  linked-diagram 
expansion. 

As  long  as  *I»0  is  the  SCF  result,  the  second-  and  third-order 
terms  in  eq.  (2)  are  determined  solely  by  C'l  double  excitation 
and  there  is  no  difference  between  RSPT  and  MBPT.  In  the 
fourth  order,  in  addition  to  double  excitations,  there  are  single, 
triple  and  quadruple  excitations  contributions.  It  is  at  this  order 
that  the  difference  between  a  truncated  C'l  (or  RSPT)  and  an 
MBPT  model  can  occur.  If  P  is  limned  to  double  excitations, 
eq.  (2)  gives  the  fourth-order  approximation  to  D-CI.  which  re¬ 
tains  unlinked,  size-inextensive  terms.  To  the  contrary  ,  the  sum 
of  the  double  excitation  diagrams  shown  in  Fig.  1  is  a  different, 
size-extensive  approximation.  On  the  other  hand,  if  double  and 
quadruple  excitations  are  both  included  in  P.  then  the  fourth- 
order  energy  obtained  is  the  same  as  the  sum  of  double  and 
quadiuple  excitation  diagrams,  or  DQ-RSPT(4)  =  DQ-MBPT(4i. 

The  coupled  cluster  method  iCCM)  |S- 1 4]  may  he  viewed  as 
a  way  to  sum  certain  categories  of  MBPT  diagrams  'o  all  orders 
and  this  also  ofters  a  somewhat  different  pity  steal  insight  ihan 
Cl  or  MBPT 

For  a  cluster  operator.  T,  one  considers  the  wave! unction 
[8-14. 20 | . 

0Cc  =  cT*o  ,-'1 

with  T  separated  into  one-body,  two-body.  etc.,  cluster 
contributions. 

f  =  f,  *  f;  -  fj  *  .  .  .  (4) 

The  various  parts  of  T  are  assumed  to  be  represented  in  tiie 
occupation  number  representation  with  the  coefficient  to  be 
determined  by  the  coupled-cluster  equations. 

For  the  example  of  Cl  quadruple  excitations.  C4  d>0.  from  the 
exponential  operator  in  eq.  (3).  we  have  the  correspondence  that 

Cx  =  fa  I  2  71  -  f,f,  *  I  :i]t:  *  I  4’f;  (5  I 

This,  in  effect,  provides  a  decomposition  of  Cl  quadruple  exci¬ 
tations  into  five  perhaps  more  physically  meaningful  com¬ 
ponents.  T\  in  j  sense  corresponds  to  two  simultaneous  inter¬ 
actions  of  two  electrons,  while  A  tends  to  represent  a  true 
four-body  interaction  [ 20 ] .  Since  the  Hamiltonian  we  are 
considering  has  no  more  than  two-body  interactions,  the  T* 
contribution  would  appear  to  be  far  more  important  than  A 
for  most  closed-shell  problems.  This  is  supported  by  pertur¬ 
bation  theory,  where  it  may  be  shown  that  all  fourth-order 
quadruple  excitation  contributions  come  from  T:  ■  with  72, 
beginning  to  contribute  in  the  fifth-order  energy  [|4|.  The  te- 
maming  terms  contain  7,  which  is  identically  zero  for  Bruecxnet 
orbitals,  and  normally  found  to  be  relatively  unimpotiant  lot 
closed-shell  systems  even  for  SCF  orbitals.  As  a  consequence, 
the  firs;  reasonable  coupled  cluster  approximation  (coupled 
cluster  doubles.  CCD)  has  the  form  [  1 0  1 3 1 
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^ccd  -  eT‘  'to  >  (6) 

By  projecting  H' Pccd  onto  'he  space  of  double  excitations, 
one  obtains  the  nonlinear  equations  of  CCD.  which  have  the 
general  form  [13) 

C> 

i  J.M 

The  coetficients  a^,  di/k,  and  »,  are  simply  combinations  of 
molecular  integials  while  the  (r;]  are  to  be  determined.  The 
equations  are  independent  of  the  energy  which  is  an  indication 
of  the  underlying  dependence  of  CCD  on  the  linked-diagram 
theorem.  A  second  important  consequence  is  that  the  CCD 
method  allows  one  to  include  the  predominant  effect  of  quad¬ 
ruple.  sextuple,  etc..  Cl  excitations  within  a  tractable  compu¬ 
tational  scheme,  while  in  the  traditional  Cl  method,  a  distinc¬ 
tion  between  the  different  components  of  C« ,  C4 .  .  .  .  is  not 
possible. 

In  the  present  work,  these  nonlinear  equations  are  solved 
iteratively.  The  first  two  iterations  of  the  linear  part  of  the 
equations  provide  the  second-,  third-,  and  fourth-order  double 
excitation  diagrams  of  Fig.  1  [14],  while  the  first  iteration  of  the 
nonlinear  part  of  the  equations  (following  a  single  linear  iter¬ 
ation)  provides  the  fourth-order  quadruple  excitation  diagrams 
shown  in  Fig.  1  [14].  The  remaining  iterations  of  these 

equations  are  not  easily  related  to  an  order-by-order  pertur¬ 
bation  approach.  The  CCD  result  is  given  as  the  converged 
solution  of  this  procedure.  (Note  that  many  solutions  of  the 
CCD  equations  are.  in  principle,  possible  [39-41],  but 
assuming  normal  convergence,  the  lowest  solution  should  be 
obtained  from  this  procedure.) 

3.  Correlation  energies  for  atoms  and  molecules 

At  the  present  state  of  development  of  molecular  theory,  as 
contrasted  with  atomic  theory,  the  use  of  conventional  atomic 


orbital  basis  sets  centered  at  the  various  atoms  or  other  speci¬ 
fied  locations  in  a  molecule,  remains  a  necessity  Some  progress 
in  numerical  methods  ior  molecuics  is  occurring  [42 1  but  these 
developments  are  likely  to  be  too  slow  to  eliminate  the  neeu 
for  basis  sets  fin  quite  some  time  As  such,  it  is  pertinent  to 
develop  information  about  what  might  be  expected  from 
MBPT  CCM  calculations  with  a  size  of  basis  set  that  can  reason¬ 
ably  be  used  for  a  number  of  chemically  interesting  problems. 
The  use  of  basis  sets  also  means  that  the  test  possible  answer  :s 
given  by  the  "full”  Cl  result  whose  agreement  with  experiment 
is  a  function  ot  the  calibre  of  the  basis  set. 

in  Tables  I  and  II  results  are  shown  for  some  atoms  and 
molecules  at  the  level  of  a  good  but  standard  contracted 
Gaussian  type  orbital  (CGTO)  basis  set.  For  C,  N.  and  0  the 
basis  consists  of  Dunning's  5s3p  contraction  [43]  of  Huzinaga  s 
primitive  9j5p  set  [44]  .augmented  by  aJ-polarization  function. 
For  H,  the  functions  consist  of  Dunning's  3s  contractions  with 
a  p-polarization  function.  The  J  and  p  exponents  are  given  else¬ 
where  [29] .  From  Table  III,  it  is  found  that  this(5s3plJ'3slpi 
basis  is  capable  of  providing  about  oO-bOfc  of  the  total  corre¬ 
lation  energy,  and  70-SO'c  of  the  valence  shell  correlation 
energy  (i.e..  neglecting  the  A'-shell  electrons  on  C.  N.  and  0). 
In  actual  values  this  amounts  to  an  error  from  about  0.02 
Hartree  up  to  as  much  as  0.20Hartree.  The  error  due  to  the 
basis  set  at  the  SCF  level  is  about  0.02  Hartree  for  H:0  and  CO. 
and  0.04  Hartree  for  CO*.  However,  in  chemistry,  all  problems 
involve  energy  differences,  and.  in  general,  much  of  the  basis  set 
error  will  cancel  to  enable  more  reliable  predictions  than  might 
be  expected  based  upon  the  accuracy  of  the  absolute  energies. 
(See  Sections  4  and  5). 

The  perturbation  energies  listed  in  Tables  I  and  II  all  refer  to 
the  standard  Moeller-Plesset  (MP)  splitting  of  the  Hamiltonian 
with  SCF  (or  Vs)  [6]  orbitals  being  used  for  the  occupied  and 
excited  one-particle  states.  It  is  apparent  that  the  second-order 
energy  provides  the  predominant  correlation  correction.  In 


Table  I.  Energies  computed  by  SCF(UHF),  fourth-order  MBPT  for  single,  double,  and  quadruple  diagrams,  and  the  coupled  cluster- 
doubles  approximation.  Basis  sets  are  (531 ).  (All  energies  in  Hartree  a.u.) 


Atoms 

fSCF 

E, 

Ef 

E? 

E? 

esdq 

DQ-MBPTI4) 

CCD 

C  CP) 

(2 S*  1) 

-  37.689  13 
3.004  2 

-0.074  50 
3.002  2 

-0.015  40 

“0.000  20 
_a 

-0.004  67 
3.001  9 

*  0.001  13 

-0.003  74 

-0.093  45 

-0.094  81 
3.0010 

N(*5) 

OS  a  1) 

-54.398  27 
4.002  7 

-0.097  54 
4.0016 

-0.01396 

-0.000  14 
_a 

-0.002  81 
4.001  1 

-  0.001  19 

-0.001  76 

-0.113  25 

-0.113  19 
4.000* 

Of‘P) 

(25  -  1) 

-74.806  98 
3.003  9 

-0.135  52 
3.002  1 

-0.009  65 

-  0.000  38 

_a 

-0.002  28 
3.001  5 

-  0.001  23 

-0.001  43 

-0.146  23 

-0.14994 
3.000  8 

0  Contribution  of  single  excitation  diagrams  to  <2S  -*■  1)  is  not  included. 


Table  il.  Energies  computed  by  SCF,  fourth-order  MBPT  for  single,  double,  and  quadruple  excitation  diagrams,  and  the  coupled 
cluster-doubles  approximation.  Basis  sets  are  (531/31).  (All  energies  in  Hartree  a.u. ) 


Molecule 

E-. 

E, 

Ef 

E? 

4? 

fSDQ 

DO-MBPTI4) 

CCD 

H.O° 

-  76.047  84 

-0.241  11 

-0.003  69 

-0.00148 

-0.003  88 

-  0.002  16 

-0.003  19 

-0.246  51 

-  0.246  66 

NH,  4 

-  56.209  34 

-0.22010 

-0.011  61 

-0.001  19 

-0.004  25 

*  0 .002  85 

-0.002  59 

-0.233  1  1 

-  0.233  40 

CH,C 

-40.206  59 

-0.189  80 

-0.01969 

-0.00097 

-0.004  89 

*  0.002  94 

-0.002  92 

-0.21143 

-0.211  ** 

CO.  d 

-187.685  91 

-0.52001 

*  0.020  53 

-0.011  77 

-0.013  36 

-  0.008  86 

-0.01627 

-0.503  98 

_  o.5oa  01 

CO  * 

-  1  12.767  33 

-0.309  82 

♦  0.004  92 

-0.006  92 

-0.009  98 

*•  0.004  90 

-0.01200 

-0.309  97 

-  0  309  t>9 

a  f^OH  =  1  308b.  9  =  1 0 2 . 4 0  > .  All  electrons  correlated. 

4  =  1.92b.  9  =  107. 3°).  All  electrons  correlated. 

c  iPch  ~  2.065b.  6  =  104.28°).  All  electrons  correlated. 
d  <^co  =  2.207b).  A'-shell  electrons  are  not  correlated. 

*  tRco  =  2.132b).  A-shell  electrons  are  not  correlated. 


r 
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Table  III  Correlation  energies  determined  by  MHPT  and  CCD  methods  Percentage  of  correlation  energy  obtained  is  given  in 
parenthesis 


All  electrons 

Valence  electrons*1 

SDO  MBPTI4) 

CCD 

1  XT6 

SDO  MBPI  (4 ) 

CCD 

1  XP6  •' 

t 

-0.0936  163  1 

-0.0948  164'") 

-  0  149 

O  0761  t?97| 

0-U773  (SI  J 

-  0  096 

N 

-  0.1  134  163') 

-  0.1 132  >63"  ) 

0.181 

O  1)949  Ity.  1 

0.0947  (75  "  i 

-  0  127 

O 

-  0.1466  l6tr  1 

-  0  1499  (61'7) 

-  0.245 

-  0.1274  168  1 1 

0  1307  1 70'"  1 

li  186 

ch4 

-0.2124  |7J7| 

—  0.2118  172":) 

-0  293 

-  0.1929  (81  ’ ) 

-  n  1933 (817  ) 

-  0.240 

Ml, 

-  0.2343  *70-7 ) 

-  0.2334  1697) 

-0.334 

-  0.2170  1 78' : ) 

0.2161  (77"  i 

-  o  280 

H.O 

-0.2480  <68'7 ) 

-  0.2467  ( 6 8"T > 

-0.365 

-  0.2297  1 757  1 

-  0  2284  1 75  ") 

■  •396 

CO 

_ 

_ 

-0.530 

-  0  3169  1767) 

-  0.3097  (747 ) 

-  o  418 

CO. 

- 

- 

-0.863 

-0.5158  1757) 

-  0.5040  (73'" ) 

-  0  692 

J  The  valence  correlation  energy  is  computed  for  CO  and  CO.  by  freezing  the  core  electrons.  Tor  other  cases,  the  second-order  pair  correlations 
involving  the  ls;  pair  are  subtracted  from  the  net  correlation  energy  given  by  the  method,  to  obtain  a  valence  correlation  estimate 
b  Reference  1 59 1.  Also  JAN  AT  Thermochemical  Tables.  National  Bureau  of  Standards.  All  correlation  energies  are  corrected  for  zero-point  '  ibrati.m 
and  relativistic  corrections. 

c  Valence  shell  correlation  energies  are  obtained  from  the  experimental  values  by  subtracting  the  I  s'- .  KK  and  KL  shell  contributions  to  the  correlation 
obtained  by  R.  k.  Nesbet.  Phys.  Rev.  175.  2  <  1968).  for  the  He  isoelectronic  sequence. 


some  cases,  such  as  CO  urul  C0;  where  multiple  bonds  are  in¬ 
volved.  second-order  can  even  somewhat  over-estimate  the 
infinite-order  CCD  energy.  Note,  however,  that  with  the  MP 
splitting,  second-order  is  usually  surprisingly  close  to  the  CCD 
result,  and  as  a  consequence,  is  probably  a  fairly  reasonable 
first  estimate  of  the  basis  set  limit. 

Alternatively,  if  the  excited-orbitals  are  determined  within  a 

‘ 1  potential  [t>]  (i.e..  a  unity  transfromation  of  the  ex¬ 
cited  orbital  space)  [24.  30.  45. 46) .  or  denominator  shifts  (6) 
are  used  to  define  a  different  splitting  of  the  Hamiltonian  (e.g.. 
the  Epstein- Nesbet  partitioning  (47|).  much  lower  second-order 
energies  are  obtained  [15].  In  many  such  cases  these  energies 
are  close  to  the  •“experimental”  correlation  energy,  but  since  the 
basis  set  limit  is  the  only  valid  objective  of  the  calculation,  and 
since  a  different  separation  of  the  Hamiltonian  or  a  unitary 
transformation  among  the  excited  orbitals  cannot  change  this 
ultimate  result,  these  modifications  will  also  make  the  third- 
and  higher-order  energies  larger  to  offset  the  large  overestimate 
in  second-order.  Since  for  most  finite  basis  sets  the  unmodified 
MP  perturbation  series  seems  fairly  well-behaved  it  seems  to  be 
preferable.  Once  CCD  or  some  similar  infinite-order  result  is 
obtained,  there  is.  of  course,  no  difference  and.  in  fact,  this 
invariance  is  essentially  achieved  in  low-order  [15.  30] . 

There  is  some  question  whether  the  perturbation  series 
would  have  the  same  observed  behaviour  if  one  actually  had  the 
"true”  U'v  and  I/N  ‘ 1  excited  orbitals  as  opposed  to  their  pro¬ 
jection  onto  a  limited  space.  This,  of  course,  is  closer  to  the 
situation  where  numerical  solutions  are  obtained  as  is  the  case 
in  atomic  calculations.  In  such  calculations,  the  Uv"'  potential 
is  also  a  useful  aid  in  performing  the  numerical  integrations 
since  the  excited  orbitals  are  constrained  from  being  as  diffuse 
as  Vs  orbitals  [6] . 

For  the  open-shell  atoms  in  Table  I  a  L'HF  reference  function 
is  employed.  Consequently  these  solutions  are  not  an  eigen¬ 
function  of  spin.  The  multiplicity  computed  from  the  transition 
state  formula.  (‘I>0iS:  i'P).  where  'P  is  the  appropriate  correlated 
wavefunction  are  also  given.  In  each  case  the  UHF  function  is  a 
good  approximation  to  the  spin  state  with  additional  improve¬ 
ment  obtained  from  higher-order  perturbation  corrections.  By 
monitoring  the  multiplicity  in  this  manner  for  open-shell  cases, 
one  can  be  confident  that  the  errors  due  to  spin  contamination 
are  not  too  great.  In  some  cases,  particularly  for  singlet  states 


where  a  UHF  function  is  employed  to  achieve  correct  separation 
at  large  internuclear  distances  to  open-shell  fragments  (e.g..  N; 
as  discussed  in  Section  5)  spin  contamination  is  sometimes 
observed  to  be  quite  large  suggesting  that  a  UHF  based  approach 
is  inappropriate. 

In  fourth-ouder.  in  addition  to  the  contributions  from  single, 
double,  and  quadruple  excitation  diagrams  given  in  Tables  I  and 
II.  there  are  the  lb  antisv mmetti/ed  triple  excitation  diagrams 
shown  in  Fig.  1.  Unlike  all  other  diagrams  in  Fig.  1  and  the 
infinite-order  CCD  model,  whose  dependence  on  the  number  of 
basis  functions  is  <  .Vs .  these  diagrams  have  an  .V'  dependence. 
This  makes  it  impractical  to  include  these  terms  without  ad¬ 
ditional  simplifications.  Judging  from  the  fact  that  the  individual 
components  of  F'?DQ  have  about  the  same  order  of  magnitude, 
one  expects  that  the  triple  excitation  terms  would  also  be  about 
the  same  size.  Furthermore,  these  terms  are  negative.  Hence,  the 
fourth-order  contribution  shown  is  actually  an  upper  bound  to 
the  full  fourth-order  result.  One  hopes,  however,  that  the 
fourth-order  triple-excitation  terms  will  have  j  relatively  small 
effect  on  energy  differences  and  may  be  safely  neglected.  This 
is  not  true  for  some  cases  [4S ]  .  however,  and  should  be 
investigated. 

Even  with  the  triple  excitation  diagrams  neglected.  /f?DQ 
usually  is  about  the  same  magnitude  as  F  3.  and  can  be  larger. 
This  is  partially  due  to  the  inclusion  of  new  ty  pes  of  excitation 
diagrams,  but  even  the  fourth-order  double  excitation  part. 
F?.  may  frequently  have  a  larger  magnitude  titan  F'3<=F'?). 
This  would  seem  to  imply  an  asymptotic  behaviour  to  the 
perturbation  series  that  could  he  a  problem.  However,  the 
CCD  model  tends  to  support  the  validity  of  the  fourth-order 
approximation.  In  CCD  one  is  summing  all  double  excitation 
diagrams  to  all  orders  and  all  the  linked  quadruple  excitation 
energy  diagrams  that  arise  from  the  disconnected  wave  function 
component  1  •  2 7"?  I *l>o  >  plus  their  mutual  coupling  [14],  Since 
this  is  done  in  an  iterative  fashion,  convergence  to  the  CCD 
result  is  observed.  From  Tables  1  and  II  it  is  further  apparent 
that  DQ-MBPT14)  and  CCD  arc  usually  quite  close. 

The  near  coincidence  of  DQ-MBPT14)  and  CCD  is  funher 
illustrated  in  Fig.  2,  where  a  plot  of  the  differences  between  the 
two  approximations  is  shown  for  a  wide  variety  of  molecules 
at  their  equilibrium  geometry.  Since  chemical  accuracy  is 
normally  thought  to  be  "  1  kcal  mole'1  (0.045 eVt,  there  is 
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Fig.  2.  Illustration  of  the  difference  in  the  energy  obtained  by  DQ- 
MBPTI4)  and  CCD  for  a  number  of  molecules.  In  the  cases  where  a  mol¬ 
ecule  is  listed  twice,  a  different  basts  set  is  used.  In  all  cases  the  basis  sets 
are  at  least  of  double-reta  plus  polarization  quality. 

little  inaccuracy  in  using  the  DQ-MBPT(4)  model  most  of  the 
tune.  At  geometries  where  the  theory  is  less  satisfactory,  such  as 
Nj  at  some  distance  beyond  equilibrium  where  a  single  deter¬ 
minant  reference  function  is  inadequate,  a  difference  of 
~  30  kcal  mole"1  is  observed  with  the  CCD  result  being  superior 
[49 1 .  Since  the  single  excitation  terms  in  fourth -order  are  not 
insignificant,  one  hopes  that  the  SDQ-MBPT(4)  results  will 
also  be  similarly  close  to  results  obtained  from  a  CCSD  model 
fi.e..  Wccsd  —  1 4>o  rl- 

From  the  foregoing,  it  is  apparent  that  the  quadruple  ex¬ 
citation  contributions  can  be  handled  within  MBPT/CCM.  The 
next  logical  question  is  their  effect  on  the  correlation  energy. 
The  magnitude  of  this  effect  may  be  assessed  by  considering  the 
difference  between  the  almost  coincident  CCD  or  DQ-MBPT(4) 
energies  and  the  energy  obtained  by  D-CI.  Some  of  these 
results  are  shown  in  Table  IV.  [Except  for  H;0  in  the  STO 
basis  set  where  we  have  done  the  all  doubles  Cl  calculations, 
the  VPD-CU2)  model  [14]  is  used  for  D-CI.  This  is  a  variational- 
perturbation  approximation  to  D-CI  that  is  correct  through 
fifth-order  in  the  energy.  In  all  cases  shown  here,  the  VPD- 
CIO  results  should  be  within  about  a  percent  of  the  D-CI 
correlation  energy  (for  H:0  VPD-CI(2)  gives  99.9%  of  D-CI 

Table  IV.  Effect  of  quadruple  excitations  in  Cl 


energy)  which  is  sufficient  to  make  a  qualitative  assessment  of 
the  importance  of  the  quadruple  excitations  on  the  energy.! 

The  molecules  m  Table  IV  are  arranged  by  the  number  of 
electrons  that  are  being  correlated.  From  size-exiensivitv 
arguments  based  upon  noninteracting  electron  pairs,  one 
expects  the  importance  of  the  quadruple  excitation  contri¬ 
butions  to  roughly  increase  with  the  number  of  electron  pairs. 
Other  factors,  such  as  the  extent  of  localization,  quality  of  basis 
set,  etc.,  also  affect  this  trend,  but  this  is  basically  what  is 
observed  from  the  present  examples.  The  double-zeta  basis  set 
used  for  benzene  is  not  capable  of  achieving  a  lot  of  correlation, 
but  a  surprisingly  large  effect  is  still  observed.  Overall,  the  error 
due  to  neglecting  quadruple  excitations  in  Cl  still  amounts  to 
4-20%  of  the  correlation  energy  even  for  relatively  small 
molecules.  In  Section  5.  it  will  be  demonstrated  that  even  an 
effect  of  4%  due  to  quadruple  excitations,  can  still  have  a 
definite  influence  on  the  accuracy  of  a  potential  energy  surface 
[19],  As  a  consequence,  it  is  frequently  unjustified  to  neglect 
such  higher  excitations  in  Cl.  and  either  many-body  methods  or 
multi-reference  based  Cl  schemes,  that  would  include  the  more 
important  quadruple  excitations  as  double-excitations  from 
doubly  excited  reference  determinants,  should  be  encouraged. 

In  table  S  several  N1BPT  CCD  calculations  are  shown  for 
H:0  in  different  sized  basis  sets.  It  is  apparent  that  the  principal 
effect  of  the  basis  on  the  correlation  energy  occurs  in  second 
order.  The  next  largest  effect  is  at  the  SCF  levei.  with  third  and 
higher-order  being  affected  by  the  choice  of  basis  set  only  at  the 
millihartree  (0.63  kcal  mole"1)  level.  In  fact,  after  second-order, 
there  is  little  to  choose  between  a  good  CGTO  basis  and  a  quite 
complete  STO  basis  set.  This  is  largely  due  to  the  fact  that  the 
STO  basis  in  this  case  correlates  the  core  electrons  much  better 
than  does  the  CGTO  basis,  but  the  higher-order  terms  are 
primarily  valence  correlation  contributions.  In  the  case  of  H;0 
good  convergence  is  obtained  even  for  the  smail  (4s2p  2s)  basis, 
but  for  Nj.  for  example,  convergence  is  much  better  in  a 
larger  basis  set. 


[(VPD-CH2»-DQ-MBPTf4»l 

Percent  error  in 

(Kcal/ mole) 

correlation  energy 

H.O  (5j3plJ/3slp)  CGTO 

6.2 

4.0 

H.ty  (5r4pM/3rlp>  STO 

7.6 

4.2 

B.H,6  (4r2pld/2slp)  CGTO 

CO6  1 5s3p  Id)  CGTO 

10.4 

6.2 

12.7 

6.5 

HCN6  I4s2p IJ/2r lp>  CGTO 

14.9 

3.2 

CO.6  i5r3plzf)  CGTO 

29.2 

9.0 

CH,CN6  i4r2pld/2a lp>  CGTO 

30.5 

11.0 

C.H,6  (4i2p'2i)  CGTO 

63.6 

19,6 

a  D-CI  is  used  for  comparison  in  this  example. 
6  Core  electrons  are  frozen  in  these  examples. 


Table  V.  Basis  set  effect  on  components  of  the  correlation  energy  of  H:0.  I  All  energies  in  Hartree  a.u.  i 


Basis 

fSCF 

ft 

e, 

E? 

E? 

E? 

ESD9 

DQ-MBPTi4i 

CCD 

i4i2p/2r)flCGTO 

-"6.0093 

-0.1378 

-0.0021 

-0.0009 

-0.0030 

-0.0008 

-0.0047 

-0.1437 

-  0.1440 

1 5r3pld' 2rlp)  CGTO 

-  76.04-3 

-0.2411 

-0.0037 

-0.0015 

-0.0039 

-  0.0022 

-0.0032 

-0.2465 

— 

I5s4p2d/3ilp>6  STO 

-  76.0642 

-0.2813 

-0.0032 

-0.0020 

-0.0043 

-  0.0032 

-0.0031 

-0.2861 

-  0.2362 

(Estimated  SCF  Limit 

—  76  0675  ) 

'Estimated  correlation  energy 

-  0  5 "C » 

“  Dunning  double-zeta  basis  set  (43 1.  The  results  for  the  valence  shell  correlation  energy  in  this  basis  through  fourth-order  have  been  reporteJ  in  ( 3-4 
6  Slater  type  orbital  basis  set  of  [511.  These  results  have  been  reported  in  (19). 
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4.  Dissociation  energies 

As  discussed  in  the  introduction,  one  of  the  primary  advantages 
of  many-body  methods  is  their  size-extensivity  property  This 
property  is  a  consequence  of  eliminating  all  unlinked  diagrams 
in  the  energy.  This  feature  has  the  corollary,  that  the  dissociation 
of  closed-shell  molecules  into  closed-shell  fragments  is  also 

size-consistent  in  Pople's  definition,  that  is  lim  £(AB)  = 

kab  *” 

£"(A)  +  F(B).  If  open-shell  fragments  are  involved,  then  in  the 
case  of  a  single  determinant  reference  function,  a  UHF  function 
is  usually  necessary  to  ensure  correct  separation.  This  imposes 
a  second  requirement  in  addition  to  using  the  linked-diagram 
theorem  to  obtain  size-extensive  and  size-consistent  dissociation 
to  open-shell  fragments. 

Several  thermochemistry  results  obtained  by  our  group 
[25,  28,  2l)|  are  listed  in  Table  VI.  In  all  reactions  except  the 
decomposition  of  HCO,  closed-shell  species  are  separating  into 
closed-shell  fragments.  Except  for  the  isomerization  energy  of 
HNC.  it  is  apparent  that  the  agreement  with  experiment  is 
excellent.  In  fact,  the  calculations  are  sufficiently  accurate  that 
we  question  the  experimental  value  for  HNC  -»  I1CN  rearrange¬ 
ment  which  will  be  the  subject  of  another  communication  [50] . 

The  basis  sets  used  in  these  calculations  are  double  zeta 
augmented  by  polarization  functions  (for  HCO  a  somewhat 
larger  basis  is  used.)  At  this  level  60-70  functions  are  required 


to  study  the  horane  containing  molecules  jnd  the  uieihyl- 
isocyamde  rearrangement  These,  along  with  oui  studies  ut 
molecular  hyperpolarizabilities  (55|.  are  the  largest  basis  MBPT 
calculations  which  have  been  made  As  mentioned  in  Section  2, 
even  though  good  basis  sets  still  permit  large  errors  in  absolute 
energies,  energy  differences  can  be  expected  to  be  much  betier. 
This  is  illustrated  by  the  icsults  of  Table  VI  Ol  course,  the  best 
results  would  be  expected  when  the  fewest  bonds  are  broken 
In  Table  VII  are  computed  the  dissociation  energies  for  the 
molecules  in  Table  11  obtained  for  decomposition  to  their 
atomic  constituents.  The  errors  in  this  case  vary  from  -  20  to 
—  80 kcal  mole ' 1  .  or  about  5  -20',7.  These  are  compared  with 
results  obtained  from  the  VPD-CII2)  method.  Even  though  the 
VPD-C1  method  is  not  size-extensive,  and  one  would  normally 
prefer  to  obtain  dissociation  energies  in  such  a  Cl  model  by 
doing  a  super-molecule  calculation  for  the  separated  species,  the 
differences  in  the  computed  dissociation  energies  are  not  that 
much  poorer  than  that  obtained  by  the  various  many -body 
models.  Agreement  may  be  further  improved  by  converging  to 
the  D-Cl  result. 

5.  Potential  energy  surfaces 

The  next  important  question  to  pose  about  many-body  methods 
is  their  reliability  for  potential  energy  surfaces.  Just  as  in  the 


Table  VI.  Comparison 
except  for  HCO  where 

of  thermochemistry  results  obtained  by 
a  tdsSpId  dslpi  basis  is  used/ 

■  SCF  and  MBPT  with 

experiment  /All  basis  sets  arc  /4s2pld  2slp: 

-  A/f  <  kcal  'mole ) 

Reaction 

Method  SCI 

MBPT/CCD 

Experiment 

2BH ,  —  B.H,** 

SDQ-MBPTU) 

18.5 

35.6 

36.6  •-  ld 

BH,  -  CO  -  H,BCO“ 

D-MBPT14) 

8.0 

20,5 

20.4  t  2d 

Bll,  -  Ml,  -  H.BMI,® 

D  MBPT(4) 

20.5 

30.1 

- 

HNC  -  HCN* 

SD0-MBPT14) 

9.7 

14.7 

110.3  r  IP' 

IINC  -  |H$|b 

SD0-MBPTI4) 

34.6 

-  30.9 

- 

CH.NC-CH  ,CN6 

SDQ-MBPTU) 

19.2 

22.  S 

23.7  :  0  14^ 

CH.NC-ICH^cl'’ 

SDQ  MBPT14) 

44 

-  40 

-  38.4" 

ti  -  CO  -  I1C0C 

CCD 

4.8 

13.7 

15.7  :  l.5h 

HCO  -  |HCO|c 

CCD 

12.8 

-  18.2 

- 

“  Reference  1 29 1 . 

6  Reference  ( 25).  Square  bracket  indicates  a  transition  state.  This  resuit  includes  a  4  kcal/molc  zero  point  correction  for  the  transition  state 
c  Reference  [  28|.  Square  bracket  indicates  a  transition  state.  This  result  includes  a  4. 8  kcal.  mole  zero  point  correction  for  the  transition  state 
d  Fehlner.  T.  P  and  Mappes.G.  W..J.  Phys.  Chem  73.  873  ( 1969). 

‘  Maki.  L.  lunpublished  results). 

'  Baehal-Vayjovec.  M.  H..  Collister.  J.  L.  and  Pritchard.  H.  O.  Can.  J.  Chem.  55.  2634  1 197_). 

*  Schneider.  F  W  and  Rabinovitch.  B  S  .  J .  Am.  Chem.  Soc.  65 .  1 794  ( 1 969). 
h  Warneck.  P  .  Z.  Naturforsch.  A26.  2047  ( 1971). 


Table  VII.  Predicted  dissociation  energies,  De.  from  different  methods  ( energies  are 


kcal /mole ,a 


Method 

CH, 

NH, 

H.O 

CO, 

CO 

CO.  - 

SCF 

324.71 

195.20 

151.14 

240.22 

170.19 

70.03 

MBPTI2) 

397  06 

272.1  1 

217.40 

349.70 

232.82 

1 16.89 

MBPTl  3 1 

399.76 

270.64 

213.66 

315.04 

214.01 

101.04 

D  MBPTl4) 

399  90 

271.54 

214.66 

317.63 

215.91 

101.73 

SD-VIBPTI4) 

400.38 

272.20 

215  35 

324.41 

219.89 

104.53 

SDO  -MBPT14) 

399.25 

271.16 

214.77 

321.11 

218.30 

102.82 

DQ-MBPT14) 

398.74 

270.41 

214.07 

314.31 

2 14  30 

100.01 

CCD 

398.10 

270.63 

211.83 

308.82 

210.94 

97,88 

Exp* 

419  49 

297.25 

232.37 

388.68 

258.85 

129.83 

VPD-CII 2) 

392.75 

264.84 

209.38 

288.05 

202.97 

85.08 

J  Dissociation  is  to  ground  slate  atoms  in  each  case.  C’H4  —  Cl  JP)  *  4Hi  :S>.  except  where  indicated  otherw  ise 
b  Reference  [  5 9 1  Also  JANAF  Thermochemical  Tables,  National  Bureau  of  Standards 
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Fig.  3  Percent  deviation  of  theoretical  from  experimental  values  for  the 
equilibrium  geometry  deft  scale)  and  force  constants  (right  scale)  of  the 
water  molecule.  Experimental  values  are  from  Hoy.  Mills,  and  Strey  (Mol. 
Phys.  24.  1265  <  1972))  and.  where  available,  from  the  revised  fit  by  Hoy 
and  Bunker  iJ  Mol.  Spectrosc.  74.  I  il979)v  Theoretical  methods 
employed  are  Hartree-Fock  'SCF>.  all  single  and  double  excitation  con¬ 
figuration  interaction  t  SD-CI).  and  the  many-body  perturbation  theory 
•  MBPT)  methods  at  third  order  (MBPT<3)).  at  infinite  order  including 
double  excitation  type  diagrams  only  (D-MBPT*  —  )).  and  coupled  cluster 
doubles  (CCD)  Technical  details  of  the  methodologies  are  given  in  (14). 
Agreement  between  theory  and  experiment  is  best  for  points  located 
nearest  the  central  horizontal  line. 

case  of  dissociation  energies,  relative  energy  differences  are 
important  and  one  hopes  that  calculations  in  sufficiently 
flexible  basis  sets  will  benefit  from  substantial  cancellation  of 
the  innate  basis  set  error.  However,  unlike  dissociation  energies 
where  only  a  few  points  are  required,  the  detads  of  the  shape 
and  curvature  of  the  surface  still  place  stringent  requirements  on 
the  basis  set  due  to  the  fine  energy  differences  that  occur  among 
the  many  computed  geometries. 

In  our  work,  we  have  studied  two  triatomic  energy  surfaces 
[19,28]  in  depth  and  various  diatomic  potential  energy  curves 
[14,  26. 27,49) .  Of  these,  our  study  of  the  H:0  molecule  [19] 
in  a  very  good  Slater  orbital  basis  set  (5s4p2d'jslp),  defined  by 
Rosenberg  and  Shavitt  [51).  is  particularly  informative.  The 
objective  of  this  study  is  to  investigate  the  quartic  force  field 
of  H:0  with  CCD.  MBPT.  SD-CI,  and  SCF  by  performing 
calculations  at  36  points,  involving  symmetric  and  asymmetric 
geometries  [19,52].  Ail  points  are  located  within  die  zero 
vibrational  displacement. 

The  results  of  the  prediction  of  the  bond  length  and  bond 
angle  and  several  of  the  quartic  constants  for  H:0  are  illustrated 
in  Fig.  3  [19[.  The  center  line  represents  the  experimental 
values  as  determined  by  Hoy,  Mills,  and  Strey  [53] ,  and  as 
revised  by  Hoy  and  Bunker  [54],  It  is  apparent  that  the  trend 
in  the  errors  is  typically  SCF  >  SD-CI  >  MBPT(3)  >  CCD  > 
D-MBPTf00).  Had  the  SDQ-MBPT(4)  results  been  included  on 
the  figure  they  would  typically  fall  between  CCD  and 
D-MBPT(~)[19|. 

The  determination  of  force  constants  from  analysis  of  the 
infra-red  spectrum  and  the  subsequent  normal  coordinate 
analysis  does  not  always  provide  an  unambiguous  force  field. 
Hence,  better  agreement  with  experiment  does  not  necessarily 
imply  the  superiority  of  the  many-body  methods  to  SD-CI.  (The 
most  appropriate  comparison  would  be  the  full  Cl  in  the  speci¬ 
fied  basis  set.)  Furthermore,  in  rigor,  CCD  should  be  superior  to 
D-MBPTl®),  but  the  neglect  of  single-  and  triple-excitation  dia¬ 
grams  in  the  CCD  model,  which  would  lower  the  energy  com¬ 
pared  to  CCD,  would  result  in  energies  much  closer  to  those 
given  by  D-MBPT(°°|.  This  is  also  supported  by  the  observation 
that  SDQ-MBPT14)  which  includes  the  single  excitation  effects, 
is  usually  ui  a  little  better  agreement  with  experiment  than  CCD 


Table  VIII.  Comparison  of  predicted  geometries  .  >}  the  form i  i 
radical  t  HCO)  with  experiment 


*?°.  A. 

1  ,  Degrees  J 

RHFa 

1.098 

1  188 

1 30.0 

I'HF6 

1.099 

1.157 

126.8 

D-MBPTia)* 

1  111 

1  188 

124.0 

Experiment0 

1.125 

1.1*5 

124.95 

0  Bruna.  P.  J.,  Bunker.  R.  J  and  Peyerimhoif.  S.  D..  J.  Mol.  Structure 
32,  217  11976). 

&  Reference  (28| 

c  Brown.  J.  M.  and  Ramsey.  D.  A..  Can.  J.  Phvsics.  53.  2252  1 19*5 1 

[  1 9 J  .  Hence,  the  slightly  better  agreement  with  experiment  of 
D-MBPT(°°)  tends  to  be  a  result  of  some  error  cancellation 
[19]. 

The  current  results,  though  certainly  not  definitive,  are 
strongly  suggestive  that  the  inclusion  of  some  effects  of  quad¬ 
ruple  and  higher  even  ordered  excitations  as  in  the  many-body 
methods,  CCD,  SDQ-.MBPT(4),  and  D-MBPTl00).  is  important 
in  obtaining  highly  accurate  potential  energy  surfaces  even  near 
equilibrium.  (This  is  not  an  isolated  observation  since  even 
Davidson’s  Cl  estimate  of  the  quadruple  excitation  contribution 
[55-53J  has  been  found  to  have  an  improved  effect  on  poten¬ 
tial  energy  surfaces  [51.60]  ).  To  some  degree,  this  may  be 
understood  by  recognizing  that  quadruple  excitations  like 
( lb- ):(3a, —  (4a,  ):(22>; ):  are  required  for  proper  dis¬ 
sociation  of  H;0  into  0<3P)  and  2H(:S)  [52] .  and  this  effect, 
is  apparently  felt  all  the  wav  into  the  well  of  the  H-0  surface 
[19]. 

The  decomposition  reaction  of  the  form'd  radical.  HCO  — 
H  ^  CO  has  also  been  studied  [28] .  including  a  determination 
of  the  transition  state  and  an  estimate  of  the  rate  constant  for 
this  system  [28.61].  This  differs  from  the  H:0  example  m  that 
the  formyl  radical  Is  an  open-shell  species  and  one  is  interested 
in  the  surface  all  the  way  to  dissociation.  For  both  reasons,  it  is 
useful  to  use  a  LHF  reference  function  in  this  study. 

Potentially,  L’HF  functions  have  a  number  of  problems.  Be¬ 
sides  their  failure  to  be  spin  eigenfunctions,  they  may  also  con¬ 
verge  to  different  symmetry  broken  solutions  at  different 
geometries  sometimes  making  it  difficult  to  use  these  functions 
in  surface  studies.  In  the  case  of  HCO  where  only  the  single 
bond  is  broken,  this  turns  out  not  to  be  much  of  a  problem, 
however. 

The  full  surface  is  studied  at  the  level  of  D-MBPTl4i  with 
CCD  results  being  determined  at  certain  critical  points.  As  :n 
H ; 0  [19],  only  minute  differences  between  D-MBPT(4)  and 
CCD  were  observed  in  this  example  [2S] .  (In  more  complicated 
cases  where  one  determinant  is  not  a  good  zeroth-order 
approximation.  CCD  and  D-MBPT(4)  can  differ  more 
significantly  [14]). 

We  find  that  for  the  full  extent  of  the  surface  the  D-MBPTi4) 
corrected  multiplicity  varies  from  2.000  at  infinite  separation  of 
H  and  CO  to  2.02  at  the  saddle  point,  decreasing  somewhat  as 
one  moves  toward  the  equilibrium  structure  of  the  radical.  As 
in  the  case  of  H-O.  the  geometry  of  HCO  shown  in  Table  VIII 
obtained  by  the  relatively  simple  D-MBPT(A)  model  tends  :o  be 
in  excellent  agreement  with  experiment  [28] . 

One  final  example  illustrates  the  important  deficiencies  of 
the  current  single  determinant  reference  function  MBPT  CCM 
methods.  It  is  well-known  that  an  RHF  function  will  not 
separate  correctly  at  large  R  for  the  vast  majority  of  molecules. 
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Pin  -J  t  ill-  and  Rill  V  potential  enerev  curves  compared  to  exper¬ 
iment  The  minima  ol  the  curves  are  superimposed.  Hie  I  I 1 1  curve  is 
lower  in  enerev  than  the  RHI  curve  at  all  R 


A  problem  of  (lux  type  demands  some  lype  ol  multirelerence 
function  based  many-body  approach  as  developed  by  Btandow 
1 3b |  .  Lindgten  [ 3 7 1  .  and  others  (3X|.  As  the  V  triple  bond 
begins  to  break,  more  than  one  configuration  (in  j  Cl  descrip¬ 
tion)  becomes  quite  important  hi  the  wavelunclion,  and  these 
extra  configurations  cannot  justifiably  be  treated  by  pertur¬ 
bation  theory.  Preferably,  all  the  configurations  that  are  re¬ 
quired  for  proper  dissociation  of  N;  I  —  IS )  should  be  included 
m  a  reference  (model)  space,  ol  p  dimension,  with  the  re¬ 
mainder  contributing  as  the  projection  manifold  of  perturbation 
theory. 

There  are  two  interrelated  problems  with  this  multirelerence 
approach.  Since  each  of  the  p~  elements  of  the  effective 
Hamdtoman  matrix  requires  a  computation  roughly  equal  to  a 
full  MBPT  CCM  calculation  with  a  single  relerence .  the  compu¬ 
tational  difficulty  rises  rapidly  the  more  configurations  that 
need  to  be  included  in  the  model  space  Also,  if  we  consider  a 
polyatomic  molecule  and  want  to  employ  all  the  configurations 
needed  to  obtain  correct  separation  for  all  possible  decompo¬ 
sition  channels,  the  dimension,  p.  can  become  quite  large.  We 
have  pursued  a  somewhat  different  approach  to  this  problem 
which  is  described  m  a  contribution  to  be  published  [ts3 1  . 
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Pie.  .'.  i  I'll  I  I  and  t  RHI  lD-MBPTl6>  and  l  RHI- 1  CCD  potential  enerev 
curves  for  N-.  The  minima  of  the  curves  are  superimposed  The 
D  -MBPTlfet  correlated  flip  curves  are  hieher  in  energy  than  the 
l)  MBPT161  RHP  curves  between  R  =  2. ft  Bohr  and  R  =  27  Bohr.  The 
1RHI1-CCD  result  extends  the  reliability  of  the  curve  over  the  (RHI  >- 
D-MBPIT61  approximation  to  somewhat  larger/?  values. 

and  consequently,  as  a  reference  function  for  MBPT  CCM.  this 
incorrect  behavior  must  eventually  manifest  itself.  A  UHF  wave- 
function  will  generally,  (but  not  always)  separate  correctly,  but 
from  spin  contamination  and  other  problems  associated  with 
different  broken  symmetry  solutions,  the  UHF  function  may 
exhibit  incorrect  behavior  on  its  path  toward  separation.  This, 
too  can  persist  even  with  correlation  included 

These  features  are  illustrated  in  the  N;  potential  curve  shown 
in  Figs.  4  and  5  [49|.  The  minimum  of  the  curves  arc  super¬ 
imposed  to  illustrate  the  differences  as  much  as  possible.  The 
spin  multiplicity  of  the  L'FIF  N;  solutions  is  found  to  be  —  3.5 
even  though  a  singlet  state  is  desired.  This  is  indicative  of  an 
enormous  contamination  from  the  low-lying  triplet. quintet  and 
septet  states.  From  Fig.  5.  it  is  evident  that  this  curious  be¬ 
haviour  of  the  UHF  solution  persists  into  the  correlated  calcu¬ 
lation.  Whereas  the  correlated  RHF  calculation  gives  quite  good 
agreement  with  the  experimental  curve  as  far  as  it  can  be  ex¬ 
pected  to  apply,  the  UHF  correlated  calculation  shows  almost  a 
reverse  curvature  only  becoming  reasonable  near  the  separated 
atom  limit. 


6.  Dipole  moments,  polarizabilities,  and  hy  perpolarizabilities 

To  address  the  question  of  MBPT  CC  D  studies  of  properties 
other  than  the  energy,  we  have  investigated  the  dependence  "I 
a  molecule  10  an  external  electric  field  |55)  The  expansion  of 
the  energy  in  an  eiectric  field  E.  is 

1*1  E)  =  MO)—  Hill,  -  1  I'aJfP,  -  1  5! p,, 

“  I  4 ' E , h) h,:P,  -  .  .  .  <M 

Summation  over  repeated  indices  is  assumed.  The  coefficients  in 
this  expansion  are  the  dipole  moment,  u.  the  polarizability  ten¬ 
sor.  a.  and  the  second-  and  third-order  (hyper)  polarizabilities, 
p  and  y  [fc>3 ) . 

The  importance  of  p  and  a  in  molecular  theory  is  well-known, 
while  the  pand  y  hyperpolarizabilities  are  ultimately  responsible 
for  the  nonlinear  optical  effects  of  molecules  [63 1.  The  most 
intriguing  aspect  of  the  latter  quantities,  is  that  the  nonlinear 
optical  effects  can  be  exploited  to  make  an  atomic  or  molecular 
gas  frequency  multiply  laser  radiation.  As  such,  a  know  ledge  of 
hyperpolarizabilities  (preferably,  as  a  function  of  frequency  1. 
can  contribute  to  a  number  of  potential  laser  devices  of  unusuai 
capabilities.  The  p  and  y  tensors  can  be  determined  by  exper¬ 
iments  employing  second-  and  third-harmonic  generation  tech¬ 
niques  [64],  and  by  using  the  Kerr  effect  [65].  However,  the 
experiments  are  difficult,  the  range  of  uncertainty  is  large,  and 
in  many  cases  values  for  p  and  y  determined  by  different  exper¬ 
imental  techniques  differ  markedly.  with  even  opposite  signs 
obtained  in  the  case  of  the  d  =  id-.-.-  +  d.-.v.v  ~  d.-> > '  hy  perpolar- 
izabihty  for  some  molecules  [64b] . 

A  first  principle  prediction  of  quantities  iike  hy  pet  polariz¬ 
abilities  places  extreme  demands  on  any  quantum  mechanical 
method.  Since  these  quantities,  unlike  the  energy .  depend  upon 
the  long-range  behavior  of  the  charge  density ,  great  care  must  N: 
taken  to  generate  a  basis  set  that  can  adequately  describe  this 
region.  Also,  unlike  a.  few  bounding  properties  can  be  usefully 
applied  to  the  higher  polarizabilities.  Finally .  correlation  is 
expected  to  be  crucial  to  determining  the  sensitive  differences  in 
the  charge  densities  that  are  requisite  in  a  predictive  theory  of 
higher  polarizabilities. 

Phvsicc  Scone  2 : 
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Employing  the  Hamiltonian. 

H  =  3!  Hut )-r  l'v(i)  -ni)-£)-r  y  Ml./)'1  -  X (9) 

T  i  >;  j 

//  =  //0(E)  -  r,j  -  rv(E)  =  H0( E)  +  T(E)  <  10) 

the  energy  in  the  presence  of  the  field  is  given  by 
tV(E)  =  ^q^jjt(E) 

*  I  (<t>0(E)iK(E)U£-0(E)-//0(E))-1RE)!',|<J>0(E)>; 

n-l  L 

(ID 

The  quantity.  Wcn^-iE)  is  the  coupled-Hartree-Fock  result. 
Equations  (9-11)  can  be  solved  by  conventional  perturbation 
theory  where  all  quantities  including  the  K1'’  potential  are 
expanded  in  a  Taylor's  series  in  the  external  field.  Considering 
an  E)  in  this  manner,  one  obtains  the  so-called  coupled 
perturbed  Hartree-Fock  (CPHF)  result  [66.  67).  The  additional 
correlation  effects  would  have  to  be  evaluated  by  considering 
the  classes  of  diagrams  with  one.  two.  and  more  interactions  with 
the  external  field.  Also,  if  one  uses  this  double  perturbation 
approach,  one  would  have  to  distinguish  between  the  CPHF 
diagrams  and  the  “true”  correlation  diagrams  m  the  linked- 
diagram  expansion  [67] . 

Alternatively,  finite  field  techniques  may  be  employed  where 
the  HF-SCF  equation  and  the  linked-diagram  expansion  are 
evaluated,  at  several  small  finite-field  strengths,  from  which 
numerical  differentiation  provides  the  various  components  of  the 
dipole  moment  and  polarizabilities.  In  the  latter  case,  the  sol¬ 
utions  of //0(E)  at  various  field-strengths  provides  the  coupled- 
Hartree-Fock  result,  which  is  equivalent  to  CPHF  [68).  The 
correlation  corrections  can  be  added  by  evaluating  the  linked- 
diagrams  as  shown  in  Fig.  1,  but  now  subject  to  field  dependent 
orbitals.  (Note,  since  all  finite-field  orbitals  are  still  determined 
in  the  appropriate  C'S(E)  potential,  the  associated  Hartree- 
Fock  cancellations  still  apply,  so  no  non-Hartree-Fock  diagrams 
need  to  be  evaluated.)  The  finite-field  procedure  has  the  advan¬ 
tage  that  all  the  coefficients  in  eq.  (8)  can  be  obtained  simply 
by  making  enough  finite-fieid  energy  calculations.  This  permits 
using  the  programs  and  methods  developed  for  the  usual  corre¬ 
lation  problem  which  usually  offer  a  more  sophisticated  level 
of  treatment  of  the  correlation  than  would  be  convenient  to 
develop  for  each  order  in  an  external  perturbation.  Also,  there  is 
a  natural  dichotomy  into  HF-SCF  and  correlation  effects.  The 
disadvantage  is  that  several  different  finite-field  va'ues  need  to 
be  considered  to  obtam  the  components  of  the  various  tensors, 
and.  of  course,  one  must  maintain  high  numerical  accuracy  in 
every  stage  of  the  computation  in  order  to  obtain  meaningful 
numerical  derivatives. 

In  Table  IX  some  results  for  the  HF  molecule  are  displayed. 
These  are  the  first  correlated  studies  of  hyperpolarizabilities 
which  have  been  reported  [3 5 j .  The  very  large  basis  set 
I6s5p4<//  5s3p)  was  chosen  essentially  by  following  the  prescrip¬ 
tion  advocated  by  Christiansen  and  McCullough  (69)  who  used 
this  scheme  to  obtain  contracted  Gaussian  orbital  basis  sets 
which  provided  good  agreement  with  completely  numerical 
coupled  Hartree-Fock  calculations  of  pz.  at zz.  IIZZZ,  and  yzzzz. 
With  our  basis  set.  agreement  with  the  numerical  CHF  results  in 
almost  perfect  for  pz.  and  a„,  within  3 %  for  and  127c  for 
7r.-«  [35  ] 

Correlation  at  the  SDQ-M3PTI4)  level  has  a  significant  effect 
on  a  and  a.  bringing  the  resuits  into  very  good  agreement  with 


Table  IX.  Dipole  moment,  polarizability,  hyperpolanzability  <>t 
HF  (basis  !bsop4d:5s3p  /  all  results  are  in  atomic  units/ 


CHF 

SDQ-MBPTi4i 

i  xpenmcnts 

A 

0.758 

0709 

070“° 

or  =  1  3<cixx  -  ayy  -  xzz  ) 

4.59 

5.58 

5  f;0 

->1  -a_ 

1.28 

1  21 

1 .3  2C 

J  ~  *^zzz  ~  ^zxx  * 

—  8.5 

-  19  9 

- 

380 

390 

- 

7yyzz 

80 

140 

- 

a  Muenter.  J  S.  and  Kemperer.  W  .  J.  Chem.  Phys.  52.  6033  •  ]9"0). 

Werner.  H.  J  and  Meyer,  W  .  Mol.  Physics  31.  855  i!976i  Zero-Point 

Correction  to  rz  included 
c  Muenter,  J.  S.,  J.  Chem.  Phys.  56.  5409 ' 19”) 

experiment.  In  the  case  of  3.  there  is  a  marked  change  due  to 
correlation,  more  than  doubling  the  CHF  result.  Similarly,  the 
components  of  y  are  increased  by  407  and  "57.  It  is  apparent 
that  correlation  is  absolutely  crucial  to  a  predictive  theory  of 
these  quantities. 

It  is  encouraging  that  the  SDQ-MBPT14)  model  seems  to  be 
rather  accurate  for  p  and  u  in  tms  example.  It  is  well  known 
that  the  single-excitation  contributions  to  properties  other  than 
the  energy  are  quite  important.  This  is  also  true  m  the  present 
case.  The  bulk  of  the  single -excitation  effects  is  included  :n  the 
initial  CHF  results  in  the  present  approach,  since  as  shown  by 
Caves  and  Karplus.  CHF  (or  CPHF)  sum  certain  categories  of 
single  and  double-excitation  diagrams  in  the  double  perturbation 
approach  to  all  orders  [67).  Even  so.  the  residual  effects  of  the 
single  excitation  diagrams  appearing  in  the  fourth-order  energy- 
are  still  significant  [35],  On  the  other  hand,  quadrupoie  exci¬ 
tation  diagrams  make  almost  no  contribution  to  ihese  properties. 
A  superior  model  to  SDQ-MBPT14),  is  probably  a  model  such  as 
SD-MBPTf°°)  or  even  CCSD  that  would  include  the  contri¬ 
bution  of  such  single  excitation  terms  to  ail  orders. 

7.  Conclusions 

We  believe  it  has  been  demonstrated  that  MBPT  CCM  can  'e 
usefully  employed  in  a  wide  variety  of  chemically  imeresrir  c 
problems.  Within  good  CC70  or  STO  basis  sets,  geometries  j-  d 
local  force  fields  may  be  predicted  very  accurately  in  most  c.sses. 
Some  properties  other  than  the  energy  can  aiso  be  obtained  to  a 
high  degree  of  accuracy.  Dissociation  energies  in  the  examples 
described  and  others  not  yet  published  [70).  can  bo  expc-.eu  ■ 
be  highly  reliable  if  only  a  single  bond  is  broken,  becoming 
so  when  comparable  dissociation  to  atomic  fragments  is  require  . 
In  the  latter  case,  however,  the  largest  error  is  due  to  basis  set 
defects  rather  than  any  inherent  weakness  in  the  methods. 

The  most  glaring  failure  of  MBPT  CCM  as  implemented  in 
this  paper  occurs  due  to  the  inapplicability  of  a  single  determi¬ 
nant  reference  function.  This  is  the  problem  in  the  N-  example 
cited,  and  will  remain  a  problem  for  a  number  of  potential 
curves  when  multiple  bonds  are  broken,  in  some  excited  states, 
and  for  the  general  open-shell  case.  The  most  encompassing 
solution  to  ibis  problem  lies  in  multireference  based  many -body 
methods  [36-38] .  A  new  approach  of  this  type  which  has  many 
desirable  features  is  descrit-’u  in  a  forthcoming  paper  [m2] . 
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ABSTRACT 


The  coupled  cluster  singles  and  doubles  model  (CCSD)  is  derived 
algebraically,  presenting  the  full  set  of  equations  for  a  general  reference 
function  explicitly  in  spin-orbital  form.  The  computational  implementation 
of  the  CCSD  model,  which  involves  cubic  and  quartic  terms,  is  discussed 
and  results  are  reported  and  compared  with  full  Cl  calculations  for  H2O 
and  BeH2-  We  demonstrate  that  the  CCSD  exponential  ansatz  sums  higher-order 
correlation  effects  efficiently  even  for  BeH2  near  its  transition  state 
geometry  where  quasi degeneracy  efforts  are  quite  large,  recovering  98% 
of  the  full  Cl  correlation  energy.  For  H2O,  CCSD  plus  the  fourth-order 
triple  excitation  correction  agrees  with  the  full  Cl  energy  to  0.5  kcal/mole. 
Comparisons  with  law  order  models  provide,  estimates  of  the  effect  of  the 
higher  order  terms  T]T2,  T^T 2 ,  T^,  and  T^  for  the  correlation  energy. 


I. 


INTRODUCTION 


Recently,  a  number  of  applications  of  many  body  perturbation 
theory  (MBPT)^-4^  and  coupled  cluster  methods  (CCM)^5”8^  to  the  ab  initio 

(9) 

calculation  of  the  electronic  structure  of  molecules  have  been  reported. 

These  applications  have  restricted  the  full  MBPT/CCM  model  to  a  fixed 
level  of  perturbation  (e.g.,  third  or  fourth  order  as  in  D-MBPT(3)  and 
SDQ-MBPT(4) )  or  to  including  all  orders  of  the  class  of  double  excitation 
cluster  operators  as  in  coupled  cluster  doubles  (CCD),  or  (poorer)  the 
linearized  L-CCD(=  D-MBPT(°°))  form.  The  methods  implemented  usually 
presume  the  use  of  Hartree-Fock  orbitals,  although  this  is  not  necessary. 

Here  we  report  the  derivation  and  the  computational  implementation  of  the 
full  CCSD  model.  The  method,  as  implemented,  uses  any  orthogonal  set  of 
orbitals  and  is  not  restricted  to  Hartree-Fock  orbitals.  In  particular, 
it  is  possible  to  use  symmetrically  orthogonal i zed  bond  orbitals  instead 
of  Hartree-Fock  orbitals  and  to  take  advantage  of  the  concomitant  reduction 
in  the  number  of  molecular  integrals  in  large  molecules  which  results  from 
the  more  localized  structure  of  the  bond  orbitals.  Additional  applications 
for  non-Hartree-Fock  orbitals,  such  as  optimizing  orbitals  so  that  the 
energy  becomes  stationary,  are  easily  envisioned. 

There  are  a  number  of  reasons  to  recomnend  the  development  of 
the  CCSD  method  as  a  basis  for  what  Pople  calls  a  "theoretical  model 
chemistry."^’^  These  criteria  propose  that  a  method(l)  should  be  "size- 
extensive,"  which  means  it  scales  properly  with  molecular  size: 
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(2)  applicable  to  a  wide  range  of  problems  within  a  single  framework; 

(3)  invariant  to  classes  of  unitary  transformations ;  (4)  efficient;  and, 

(5)  able  to  correctly  separate  a  molecule  into  its  fragments. 

CCSD  is  a  many-body  method  built  upon  the  linked-diagram  theorem. 

Hence,  it  is  size-extensive;  and,  in  particular,  CCSD  gives  the  correct  result 
for  the  characteristic  problem  of  separated  electron  pair  bonds  as  in  a 
lattice  of  N  noninteracting  H2  molecules.  As  long  as  a  single  determinant, 
which  need  not  be  a  restricted  or  unrestricted  Hartree-Fock  function, is  a 
reasonable  starting  point,  CCSD  is  applicable  to  most  problems  without 
modification  or  special  symmetry  conditions. 

Furthermore,  CCSD  is  invariant  to  any  transformation  among  the 
excited  orbitals  or  the  occupied  orbitals,  respectively .  CCSD  is  not 
generally  invariant  to  transformations  that  mix  occupied  and  unoccupied 
orbitals  among  themselves.  However,  for  the  special  case  of  noninteracting 
separated  electron  pairs,  it  is  even  invariant  to  such  general  transformations. 
This  follows  from  the  fact  that  CCSD  is  the  "full"  Cl  (i.e.,  all  possible 
contributing  n-tuple  excitations  of  n  electrons)  for  such  a  model  problem. 

Since  interpretations  of  chemistry  are  largely  based  upon  the  concept  of 
separated  electron-pair  bonds  this  is  a  very  desirable  aspect  of  the 
CCSD  model.  In  a  real  molecule,  different  choices  of  the  molecular 
orbitals  will  give  different  energies,  but  we  would  expect  a  smaller 
effect  due  to  such  transformations  for  CCSD  than  in  less  complete  models. 

It  will  be  interesting  to  see  if  localized  orbital  models  will  be 
approximately  invariant. 

The  condition  of  efficiency  also  recommends  CCSD  since  the 
treatment  of  the  electron  correlation  grows  no  more  rapidly  than  the  sixth 
power  of  the  number  of  basis  functions,  M®.  Thus,  CCSD  involves  no  more 

- . 
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coefficients  than  the  configuration  interaction  single  and  double  excitation 
model,  SD-CI;  and,  in  principle,  it  is  no  more  time  consuming.  This  is 
true,  even  though  CCSD  includes  contributions  of  quadruple  excitations, 
as  well  as  additional  parts  of  the  triple  and  higher  excitations  that 
arise  due  to  disconnected  products  of  single  and  double  excitations.  Any 
attempt  to  exceed  the  CCSD  mcdel  by  including  higher  categories  of  excitation 
operators  such  as  the  connected  triple  excitations,  T^,  will  invariably 
result  in  a  model  where  the  number  of  operations  would  asymptotically  rise 
more  rapidly  than  M®. 

The  condition  of  correct  separation  depends  upon  the  reference 

function  as  well  as  the  degree  of  correlation  included  and  generally 

requires  the  resolution  of  degeneracy  problems.  Full  Cl  with  a  single 

reference  function  obviously  separates  correctly  even  for  poor  choices  of 

reference  function.  For  less  complete  correlation  models,  the  relative 

importance  of  the  reference  function  and  the  correlation  corrections  is  not 

yet  determined,  often  recommending  multi  reference  techniques.  A  single  UHF 

function  will,  in  general,  separate  correctly,  though  it  can  suffer  from 

large  amounts  cf  spin-contamination  causing- an  erroneous  behavior  of  a 

potential  energy  surface.  With  any  choice  of  reference  function,  CCSD  will 

certainly  go  farther  toward  achieving  correct  separation  than  SD-CI,  so  we 

might  expect  a  higher  level  of  applicability.  For  example,  CCD  is  known 

(12  131 

to  correctly  handle  some  severely  degenerate  problems,  ’  and  in  this  paper, 
CCSD  is  similarly  shown  to  resolve  two  highly  degenerate  problems  without 
resorting  to  multireference  function  techniques.  Consequently,  CCSD  offers 
a  potentially  attractive  model  for  a  "theoretical  model  chemistry"  that  can 
often  even  achieve  correct  separation  or  resolve  some  kinds  of  degeneracies 
despite  employing  a  single  reference  configuration. 
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In  the  following  section  we  will  review  coupled-cluster  theory 
and  the  singles  and  doubles  approximation.  Section  III  sunmarizes  the 
CCSD  equations  and  Section  IV  sketches  the  method  of  implementation. 
Finally,  Section  V  compares  CCSD  calculations  at  the  double  zeta  level 
on  H2,  H^O,  and  Be!^  with  full-CI  calculations. 
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THE  SINGLE  AND  DOUBLE  EXCITATION  APPROXIMATION  IN  COUPLED  CLUSTER  THEORY 

Coupled  cluster  (CC),  many-body  perturbation  theory  (MBPT),  and 
configuration  interaction  (Cl)  are  methods  designed  to  solve  the  Schrodinger 
equation,  and  consequently,  all  can  be  interrelated.  The  primary  difference 
is  in  how  higher  order  excitation  configurations  are  handled,  with  the  con¬ 
sequence  that  many-body  methods  are  size-extensive.  In  particular,  CC  theory 
can  be  viewed  as  a  way  to  sum  certain  categories  of  many-body  perturbation 
theory  diagrams  to  all  orders.  While  such  a  viewpoint  has  advantages,  we 
will  discuss  the  CCSD  approximation  from  the  viewpoint  of  an  exponential 
representation  of  the  exact  wavefunction.  That  is,  we  will  write  the  exact 
solution  to  the  Schrodinger  equation  as  an  exponential  of  the  cluster  operator 
operating  on  a  reference  function, ^ 


^exact  ~  ^CC 


eT  , 


(1) 


where  is  a  single  determinant  and  T  is  a  cluster  operator  which  is 
usually  separated  into  one-body,  two-body,  etc.  cluster  contributions. 


T  =  Ti  +  T2  +  T3  +  ...  (2) 

The  various  parts  of  the  cluster  operators  are  represented  as  expansions 
of  second  quantized  excitation  operators  and  the  problem  of  determining  T 
is  reduced  to  the  problem  of  finding  the  expansion  coefficients  of  the 

A  A 

second  quantized  operators.  For  T-j  and  T^,  the  expansions  are 

h  ’  E  ‘1  *1  ai 
1 


(3) 
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and 


z 

i>j 


a>b 


(4) 


where  we  have  adopted  the  convention  that  the  lower  case  Roman  subscripts 
(superscripts)  i,  j,  k,  t ,  ...  (a,  b,  c,  d,  ...)  refer  to  orbitals  which 
are  occupied  (unoccupied)  in  the  reference  determinant.  The  undetermined 

3  9b  4.  a. 

coefficients  are  { t . >  and  (t. .}  while  (a, a.}  and  {a  a.  a,  a.}  are  second  auantized 
i  ij  a  i  a  I  b  j 

sets  of  single  and  double  excitation  operators,  and  t®*-  is  antisymmetric 


(i.e. 


$>• 


tab  .  .tab  =  -tba 

ij  JP  iJ 

The  chief  advantage  of  CC  over  Cl  can  be  easily  demonstrated  by 


using  the  expansion  of  an  exponential  operator, 


T 

e 


1  1 

i  +  T  +  7j-r+^-r  + 


(5) 


to  bring  the  CC  cluster  operators  into  a  formal  correspondence  with  the  Cl 
excitation  operators.  In  the  Cl  method,  the  exact  wave  function  can  be 
written  as  a  linear  combination  of  excitations  from  a  reference  determinant. 


exact 


=  fCI  -  (1  + 


+  C2  + 


V 


(6) 


where  C.  is  a  sum  of  i -fold  excitations  with  coefficients  which  must  be 
determined.  The  correspondence  between  Ci 's  and  TVs  is  established  by 
using  equations  (2)  and  (5)  in  (1)  to  produce  an  expanded  equation 


exact 


=  (1  +  T,  +  T, 


+  f3  +  -+  ?1  +  TiW1  +  ->iv 


(7) 


Terms  with  the  same  total  excitation  level  are  collected  together  and  equated 
to  the  Cl  coefficient  of  the  same  excitation  level.  For  example,  the  Cl 
quadruple  excitations.,  C^$,  correspond  to  the  following  sum 


1  "2 

=  T  +  — 

4  2  2 


+  T1T3 


1  *2* 

2  12 


1 

4! 


(8) 


Thus,  CC  can  be  regarded  as  a  way  to  decompose  Cl  coefficients  into  other, 

possibly  more  physically  meaningful  terms. In  the  case  of  quadruples, 

a  substantial  body  of  work  has  indicated  that  the  largest  component  of  C4 

comes  from -^2  in  equation  (8),^"^  From  a  computational  point  of  view, 

*2 

this  is  an  important  observation  because  the  effects  of  T2  can  be  included 
using  algorithms  where  the  work  is  proportional  only  to  the  sixth  power  of 
the  basis  set  size  instead  of  proportional  to  the  eighth  power  of  the 
basis  size. 

The  CCSD  method  is  an  approximate  CC  method  in  which  the  exact 
wave  function  (c.f.  equation  (1))  is  approximated  by  truncating  the  expansion 
of  T  (c.f.  equation  (2))  after  T2.  Thus, 


V  —  Y 

exact  -  CCSD 


(9) 


The  coefficients  which  must  be  determined  are  just  those  given  in  equations 
(3)  and  (4),  and  the  number  of  unknown  coefficients  in  the  CCSD  approximation 
equals  the  number  of  coefficients  in  SD-CI.  Thus,  the  level  of  computational 
effort  required  in  the  CCSD  model  is  expected  to  be  comparable  to  the  level 
of  effort  required  for  the  SD-CI  model .  However,  as  indicated  in  equation  (8), 
the  CCSD  approximation  incorporates  parts  of  the  Cl  quadruple  excitation 
terms,  namely,  jT^,  T^T^,  and  ^T^,  and  it  does  so  more  economically  and 


compactly  than  a  SDQ-CI  calculation  can.  There  is  evidence  based  upon 
perturbation  theory,  previous  calculations  and  physical  grounds  that  the 
single  missing  term,  T^,  is  usually  not  needed  to  achieve  accurate 
calculations.^  ^  In  addition  to  including  the  effects  of  Cl-type  Quadruple 
excitations,  the  CCSD  model  also  incorporates  some  of  the  effects  of  triple 


excitations.  Associating  C^  with  cluster  terms,  we  find 


T1T2 


+ 


1_ 

3! 


(10) 


■  -  ~3 

Thus,  the  CCSD  model  incorporates  the  T^T^  and  T^  components  of  Cl  type 
triple  excitations.  Unfortunately,  when  Hartree-Fock  orbitals  are  used. 


perturbation  theory  indicates  that  the  dominant  contribution  to  C3  usually 
comes  from  T.,  which  is  not  included  in  the  CCSD  method.  On  the  other  hand, 

if  non-Hartree-Fock  orbitals  are  used,  so  that  T-j  is  large,  then  the  dominant 
contributions  to  C3  can  come  from  T-j^  and  T^.  In  practice,  we  have  found 
that  these  terms  do  become  important  in  bond  breaking  processes  which  one 


is  usually  tempted  to  describe  with  a  multireference  approach.  They  seem 


to  remain  unimportant  for  simple  closed-shell  molecules  at  their  minimum 


energy  geometries.  As  in  the  case  of  the  C^  components  inherent  in  CCSD, 

a 

the  disconnected  C3  components  can  be  computed  with  algorithms  in  which  the 
work  is  proportional  to  the  sixth  power  of  the  size  of  the  basis  set. 
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III.  CCSD  EQUATIONS 

The  procedure  for  deriving  the  CCSD  equations  have  been 
previously  reported  in  diagrammatic  form^)  and  approximations  to  these 
equations  in  orbital  form/^  using  elegant  and  compact  notations. 

In  some  cases  the  reported  equations  have  been  restricted  to  Hartree-Fock 
orbitals  or  approximations  were  made  in  the  equations.  In  this  section 
we  present  the  complete  set  of  equations  satisfied  by  the  coefficients 
which  define  T-|  and  T ^  in  the  CCSD  method.  These  equations  are  applicable 
to  any  set  of  orthonormal  spin-orbitals.  In  particular,  the  eauations  are 
applicable  to  RHF  and  UHF  as  well  as  non-Hartree-Fock  reference  determinants. 
The  equations  are  derived  algebraicly  for  this  work  using  the  conventional 
procedure  beginning  with  the  Schrodinger  equation 

/ n  r  )  w  -  E  )e  ^  ^1$  >  =  0  (11) 

(H  “  ECCSD  -CCSD  uCSD  e  1  O 

and  projecting  onto  a  set  of  functions,  <®Q| >  and  {<ijl}’  such  that 

a  set  of  equations  sufficient  for  determining  the  ta  and  t*!j  coefficients 

results. 


«>0I(H  -  ECCSD)|eV\>  =  ° 


_T1+T2 


<  -|(h  -  ECcsD)|e  V  =  0  for  a11  i’  3 


(12) 

(13) 


.  t,+t7 

-  ECCSD)|e  =  0  for  a11  a>’  b 

Next,  the  exponential  is  expanded  using  equation  (5)  and,  using  the  fact 
that  H  contains  no  more  than  two-electron  operators,  we  have 
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W2. 


ECCSD,1+T1+T2'TT1  '  V 


(12a) 


a 
< . 


if-  *° 


(13a) 


-  eccsD^+T1+T24T?+TiT24T!4TfT241TH1Tt>  i  =0=-* 0  o*» 

Evaluating  (12  )  in  terms  of  the  amplitudes  ta  and  ta*?,  then  gives 


eccsd 


<$  I H [  <S>  >  +  £ 

o  o  ,• . 


f . 
ia 


ta  +  £  <i j | |ab> 
1  1>j 

a>b 


(taD  +  tat°  -  tatD)  (15) 
ij  i  J  l  J 


where  f]a  =  <f  1 H 1 5Q>  and  <i  j  J  |ab>=  j  ;<i  (1  )Xj  (2)  —  -  -  ;<a(1  )xb(2)dxT  dx2 
Finally,  equation  (15)  is  substituted  into  equations  (13)  and  (14) 
causing  terms,  which  would  make  unlinked  contributions  to  the  energy,  to 

cancel  and  eliminating  the  explicit  E  dependence  of  the  t  coefficients. 

The  other  important  step  in  the  derivation  is  to  observe  that  the  T-j 

equation  may  be  factored  from  the  equation  for  the  T,,  amplitude, 
as  discussed  more  fully,  below.  Evaluating  the  resulting  second-quantized 
matrix  elements  with  KOMMUTE*;^  a  computer  program  for  determining 
matrix  elements  between  determinants,  and  carrying  out  the  above 
simplification  results  in  the  two  equations  satisfied  by  the  amplitudes 

A  A. 

in  T-j  and  These  are  presented  in  Tables  (1,  2,  3).  After  solving 
the  equation  of  Tables  1-3  the  energy  is  given  by  Eq.  (15). 

An  important  aspect  of  the  coupled  cluster  equations,  in  general, 
and  the  CCSD  equations  in  particular  are  their  fully  "connected"  diagrammatic 
form.  All  unlinked  diagrams  (diagrams  that  contain  a  closed  disconnected 
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part)  cancel  with  the  energy  in  Eqs.  (13-14),  thereby  ensuring  the  size- 

extensivity  of  the  model/  '  Disconnected  diagrams  (diagrams  that  have 

open  disconnected  parts)  would  still  remain  in  the  general  case.  However, 
(18) 

Lindgren'  '  has  proven  that  in  coupled  cluster  theory,  only  connected 
diagrams  need  be  considered  for  the  various  Tn  amplitudes.  This  feature 
becomes  transparent  in  the  direct  algebraic  derivation  presented  here,  since 
all  disconnected  diagrams  Stained  in  the  straightforward  evaluation  of 
Eq.  (14)  correspond  to  a  single  particle  amplitude,  t*  multiplied  by  the 
equation  of  Table  1,  which,  of  course,  vanishes.  Diagrarmatical ly ,  these 
disconnected  terms  are  of  the  form 


x 


a  T,  amplitude  is  signified 


hv 


> 


(19) 

A  recent  thesis  by  Cullen.  presents  the  diagrammatic  deviation 
of  the  CCSD  equations  with  a  thorough  enumeration  of  the  diagrams. 
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Several  points  about  the  expressions  for  the  T-j  and  T2  amplitudes 

should  be  made.  First,  at  this  point  we  have  not  yet  prescribed  a  method 

for  solving  these  equations.  Second,  the  equations  are  nonlinear  and  coupled, 

a  featur ?  they  share  with  the  Hartree-Fock  or  MCSCF  equations.  Third,  the 

equations  are  a',-'*'ti c  in  T-j  but  only  quadratic  in  T_.  The  quadratic  non- 

linearities  are  similar  to  the  nonlinearity  that  is  implicit  in  the  Cl 

secular  problem  (see  Appendix).  Although  the  quartic  nonlinearity  can  be 

a  problem,  anyone  who  has  solved  Hartree-Fock  on  MCSCF  problems  is  unlikely 

to  be  deterred  by  the  comparatively  low  level  of  nonlinearity  in  the  CCSD 

equations.  A  fourth  feature  is  that  the  computational  effort  required  to 

evaluate  each  amplitude  grows  asymptotically  no  more  rapidly  than  the  sixth 

power  of  the  number  of  basis  functions.  In  a  system  with  n  occupied  orbitals 

and  N  unoccupied  orbitals  the  CCSD  computation  time  for  very  large  basis 

2  4 

sets  and  many  electrons  will  grow  only  as  n  N  ,  wh^'ch  is  the  same  as  for 

SD-CI.  Enhancements  to  CCSD  incorporating  all  excitations  in  a  class  of 

higher  excitations  (e.g.  Tj)  will  result  in  algorithms  in  which  the  time 

required  grows  as  tne  seventh  power  of  the  size  of  the  basisl20^ 

The  final  point  to  be  made  concerns  the  form  of  the  equations  for 

the  T2  amplitudes  in  Tables  2  and  3.  These  equations  have  been  written  in  a 

A  1  "2 

form  which  emphasizes  the  similarity  between  T2  and  j  T-j  terms.  Thus,  it 

-  2 

is  possible  to  implement  terms  containing  T^  using  subroutines  already 
written  for  T2  terms. 
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Implementation  of  the  CCSD  Method 


At  the  time  when  we  reported  implementing  the  CCD  model  we 
briefly  outlined  the  procedures  we  used  for  solving  the  nonlinear  equations. 
We  omitted  many  details  which  we  felt  were  self-evident;  however,  our 
subsequent  experience  has  convinced  us  that  a  more  detailed  description  of 
the  algorithms  should  be  disclosed. 

First,  we  wish  to  define  what  we  do  not  do.  If  we 
collect  the  T-j  and  1^  coefficients  into  a  vector,  X,  and  define  the  arrays 
A,  B,  C,  D,  and  E,  it  is  possible  to  write  the  equations  in  a  qeneral 
tensor  form 


Ai +  ]  BUXJ +  ,EJ  cfjkV  +  lit  DmtYA +  * 0 


In  principle,  equation  (is)  can  be  solved  using  standard  approaches  after 

constructing  the  matrices  A,  B,  C,  D,  and  E.  Although  we  find  equation  (16) 

occasionally  useful  in  discussing  properties  of  the  coupled  cluster 
(131 

equations;  the  large  dimension  of  the  B,  C,  D,  and  E  arrays  make 
implementation  of  practically  useful  algorithms  predicated  upon  the 
construction  and  manipulation  of  these  arrays  impossible.  Consequently, 
the  programs  never  construct  these  arrays;  however,  it  is  possible  to 
think  of  our  programs  as  producing  the  result  of  B,  C,  D,  and  E  operating 
on  X  without  explicity  constructing  the  arrays. 

The  methods  used  within  the  programs  can  be  most  simply  explained 
by  referring  to  the  eauation  in  Table  1  as  an  example.  The  first  step  is  to 
rearrange  the  equation  into  an  explicit  equation  for  the  coefficient  t^. 
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Equation  (19)  illustrates  one  additional  feature  of  the  method  we  use  to 
solve  the  CCSO  equations.  Namely,  we  usually  choose  to  simultaneously  iterate 


< 
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the  equations  for  and  T ^  coefficients  so  that  at  the  end  of  the  nth 
cycle  we  have  all  tcl(n)  and  tf^n).  Of  course,  it  is  possible  to  iterate 
separately  for  the  T-j  coefficients  with  frozen  coefficients  and  then  iterate 
for  T^  coefficients  with  fixed  T i  •  We  usually  find  that  T-j  and 
coefficients  are  sufficiently  coupled  that  such  an  iteration  method  is 


uneconomical . 

Another  point  which  equation  (19)  illustrates  is  that  we  retain 

the  nonlinearity  of  the  CC  equations  throughout  their  solution.  Usually, 

we  do  not  choose  to  solve  a  linear  approximation  before  introducing  the 

nonlinear  terms.  Nor  do  we  choose  to  use  a  Newton-Raphson  method  to 

achieve  rapid  convergence  because  the  N-R  method  requires  evaluation  of  a 

gradient  matrix.  The  size  of  the  gradient  matrix  would  be  too  large  to 

handle  conveniently.  Instead,  we  use  a  reduced  linear  eouation  method  to 

(131 

accelerate  convergence.  ' 

In  writing  down  equation  (17),  we  choose  to  use  the  terms 

containing  diagonal  Fock  matrix  elements  to  solve  for  t^  .  This  choice 

implies  that  the  first  few  iterations  of  equation  (17)  beginning  with  T-j  =  0 

correspond  to  a  perturbation  solution  of  the  CC  equations  using  a  Mtfller- 

Plesset  partitioning  of  the  Hamiltonian,  which  has  been  shown  to  normally 

(22) 

have  better  convergence.  If  a  linear  approximation  to  equation  (17)  is 
made,  then  all  iterations  can  be  made  just  as  in  perturbation  theory. 

Although  the  Moller-Plesset  partitioning  has  been  shown  to  be  superior 
for  perturbation  methods  based  upon  Hartree-Fock  orbitals,  there  are 
times  when  it  is  clearly  inappropriate  such  as  when  RHF  orbitals  are 
used  for  an  open-shell  configuration  and  f^  =  f^.  In  that  situation, 

(f ^  -  fdcJ)  ^  is  indeterminant  and  we  have  a  so-called  "dangerous  denominator." 
Fortunately,  in  this  case  the  dangerous  denominator  problem  is  artificial 

and  can  be  eliminated  by  electing  to  add  an  arbitrary  constant  times  to  both 

r 


v 

j 
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sides  of  the  equation  (Table  1)  and  then  solve  for  t£.  A  simple 
rearrangement  to  other  forms, like  the  Epstein-Nesbet  partitioning  for 
example,  will  not  usually  alleviate  the  denominator  singularity, 
any  circumvention  does  not  resolve  the  denominator  singularity.  Of  course, 
mul tidetermi natal  description  implied  by  the  dangerous  denominator. 

Throughout  this  paper,  the  equations  have  been  written  in  spin 
orbital  form.  However,  prior  to  implementation  in  the  computer  program, 
the  equations  are  rewritten  and  the  spin  factors  are  specifically  included. 
For  example, 

a  b 


Z  <id|jab>t*5 

a>b  u 

i 


z  <i  d  | [a  b  >  t.a„“ 
,  h  a  a  a  a  i  l 
a  >D  a  a 

a  a 

ia 


+  1  ‘'b'.iv.’  'u  120 > 

a.'-ba  8  “ 

\ 

where  i  refers  to  the  spatial  function  of  the  ith  spin-orbital  ard  b  orbitals 

are  numbered  higher  than  a  orbitals.  If  na  (n  )  is  the  number  of  occupied  orbitals 

with  a(s)  spin  and  N^Njis  the  number  of  unoccupied  orbitals  with  a(s)  spin 

2 

then  the  factoring  by  spin  reduces  the  sum  on  the  left  side  from  (n  +  n  ) 

a  3 

^  op  2 

(N  +  N  ) J  to  n  N  (N+U/2 +n  n„N  N  operations.  The  equations  are  also 
a  6  aa'a'/'agagr  n 

analyzed  to  reveal  simplifications  which  result  when  orbitals  are  spin 
restricted  so  that  a  and  8  components  have  the  same  space  functions.  Thus, 
in  the  spin-restricted  problem  tda  is  evaluated,  but  t^j®  (which  equals  t^a) 

la  la 

is  not  evaluated.  Also  terms  like  the  second  term  on  the  right 
hand  side  of  equation  (20)  simplify  as  follows  in  the  RHF  case. 
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z 

a  b 

3  a 


a  b 

<i  d  |ab  >  t.6»a 
8  a 1  6  a  1  l 
3  a 


2  .  \  <f  a«|a  b  ”  *1  £  (1-  sab/2) 

a-j  >D  a 
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2 

Consequently,  the  number  of  operations  required  drops  from  nangNaNg  to 

p  p 

n  N  (N  +  l)/2.  In  all  cases  we  have  been  able  to  implement  the  spin- 

ex  a  a 

restricted  sums  with  the  same  subprograms  as  used  for  the  implementation 
of  the  spin-unrestricted  sums  by  changing  loop  limits  and  inserting 
appropriate  factors  of  two.  As  a  result  of  explicitly  treating  the  spin, 
the  work  involved  in  evaluating  our  equations  is  essentially  the  same  as 
if  we  had  adopted  a  spin-adapted  formulation  while  retaining  the  flexibility 
of  removing  spin  restrictions. 

In  addition  to  an  explicit  treatment  of  spin  summations,  we  also 
factor  terms  containing  products  of  three  terms  into  intermediate  partial 
sums  containing  just  two  terms.  Thus, 

E  <1j||ca>  tj  t}!  »  Z  t£  a  d  (22) 

i>j  L  1J  c  4  CQ 

ca 


where  acd 


E  <1J 
1>J 
a 


|ca>  t 


da 

ij 


The  intermediate,  acd. 


is  computed  and  stored  as  partial  sum  before  completing  the  evaluation  by 
summing  over  the  T-j  terms.  As  a  result  of  this  factorization,  the  work  required 
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2  3 

to  evaluate  the  above  term  is  reduced  from  n  (n-l)N  /2  operations  to 
2  2  2 

n  (n-l)N  /2  +  nN  operations.  This  is  exactly  the  same  type  of  simplification 
required  to  reduce  the  two-electron  integral  transformation  to  an  (n  +  N)^ 

g 

process  from  an  (n  +  N)  process;  and  like  the  explicit  spin  summations,  it 
always  has  been  used  in  our  MBPT/CCD  programs. (20-24) 

Our  treatment  of  symmetry  within  the  program  occurs  at  two  levels 
and  is  intertwined  with  the  choices  we  made  for  integral  and  coefficient 
storage.  Hence,  a  brief  explanation  of  our  integral  storage  technique  is 
required.  Before  we  adopted  any  storage  scheme  we  first  proposed  some 
design  goals  for  the  program: 

(1)  The  storage  scheme  adopted  should  facilitate  writing  the 
program.  That  is,  it  should  be  easier  to  write  the  required 
subroutines  and  easier  to  ensure  the  correctness  of  the 
routines. 

(2)  Every  step  in  the  evaluation  of  the  terms  required  should 

be  fully  factored  and  no  term  should  require  more  than 

2  4  4  2  3  3 

n  N  ,  n  N  or  n  N  operations. 

(3)  Access  to  data,  both  in  memory  and  on  disk,  should  be 
sequential  in  the  inner  loops.  No  input/output  would 
ever  occur  within  the  two  innermost  loops.  Random 
access  to  records  on  disk  would  be  presumed,  but  primarily 
used  to  position  subfiles  which  would  then  be  read 
sequentially.  Because  of  the  sequential  access  through 
memory  the  program  would  be  ideally  suited  for  vitual 
memory  computers. 


•3 
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(4)  The  central  memory  requirements,  should  be  proportional 
to  the  square  of  the  number  of  basis  functions.  Thus,  if 
the  computer  has  enough  memory  to  solve  the  Hartree-Fock 
problem,  then  there  should  be  enough  memory  to  solve  the 

CC  problem.  Unlike  most  Cl  programs,  we  do  not  assume 

that  all  of  our  coefficients  could  fit  into  memory.  Instead, 
we  choose  to  hold  only  a  single  distribution  of  coefficients 
in  memory  and  to  carefully  manage  the  concommitant  increase 
in  input/output  by  working  with  fully  ordered  integrals  and 
coefficients.  Consequently,  we  were  able  to  carry  out 
frozen  core  CCD  double  zeta  calculation  on  benzene  (60  MO's) 
on  a  VAX  11-780  using  a  physical  memory  working  set  of 
256  K  bytes  while  paging  at  a  relatively  slow  rateJ11) 

(5)  Symmetry  zeroes  and  accidental  zeros  should  be  treated 
transparently  and  on  equal  footing  where  possible. 
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Our  approach  to  implementing  these  six  goals  centers  on  the 
observation  that  most  terms  in  Tables  1-3  can  be  written  as  simple  scalar 
products.  For  example,  by  combining  ab,  cd,  and  ij  into  single  labels 
c,  f,  and  k  we  can  rewrite 

E  <ab||cd>t^  (25) 

cd 

as  a  simple  matrix  product  z  Tf(< 

If  all  <ab[|cd>  and  t^  are  sorted  onto  a  random  access  device  so  that 

the  label  cd  identifies  a  record  containing  <ab||cd>  for  all  ab  and  the 

label  ij  identifies  a  record  containing  t^  for  all  ab,  then  the  summation 

*  J 

indicated  in  (24)  is  easy  to  perform,  especially  if  the  integrals  are 
also  ordered  within  each  record.  Furthermore,  the  sums  can  be  set  up  so 
that  only  integrals  of  <ab||cd>  type  and  fr  coefficients  of  t^  type 

must  be  in  memory  to  generate  a  given  t?1?  .  If  integrals  are  ordered  in 

*  J 

each  record,  exact  zeroes  and  integrals  smaller  than  a  given  threshold 
can  be  removed  and  a  skip  count  indicating  the  distance  between  labels 
can  be  packed  into  the  integrals. 

As  a  result,  to  simplify  programming, we  opt  to  sort  the 
molecular  orbital  integrals  coming  out  of  the  two-electron  transformation 
into  antisymmetrized  combinations  with  Dirac  type  labels  and  to  store 
them  in  random  access  subfiles  according  to  the  number  and  location  of 
occuppied  indices.  Within  each  subfile,  two  integral  labels  are  used  to 
specify  a  record  containing  all  integrals  of  that  type.  Zero  integrals 
and  approximately  zero  Integrals  are  not  stored  in  the  records  and  skip 
counts  are  inserted  to  keep  track  of  these  integrals.  While  these  steps 

4 

introduce  some  inefficiency  and  redundancy  at  the  N  level,  we  are  willing 
to  make  these  kinds  of  sacrifices  to  speed  up  the  N6  processes. 
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Returning  to  the  treatment  of  symmetry,  we  see  that  at  the 
lowest  level  symmetry  is  implemented  by  removing  all  zeroes  from  the 
anti  symmetrized  integral  and  coefficient  lists.  While  it  is  important 
to  remove  zeroes  for  integral  and  coefficient  storage,  the  principal 
benefit  occurs  during  the  construction  of  the  scalar  products  where  a 
routine  like  that  shown  in  Table  4  can  be  used  to  perform  a  sparse  scalar 
product.  Thus,  the  inner  most  loops  are  implemented  with  a- sparse  scalar 
product  subroutine  which  is  driven  by  the  indices  implicit  within  the 
anti  symmetrized  integral  list.  Since  indices  which  are  zero  by  symmetry 
never  appear,  the  loops  effectively  run  only  over  symmetry  indices.  To 
avoid  generating  a  term  which  is  zero  by  symmetry,  we  use  a  symmetry  template 
in  the  outermost  loops.  Thus,  the  target  arrays  contain  a  bit  flag 
which  indicates  whether  the  sum  is  zero  by  symmetry.  In  effect,  the  outermost 

loops  run  over  all  orbital  indices,  but  the  inner  loops  are  skipped  altogether 
if  the  evaluation  term  must  be  zero  by  symmetry. 

This  completes  the  discussion  of  the  computational  considerations. 
The  following  section  discusses  some  applications  of  the  CCSD  model  to 
molecules. 
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Our  first  application  of  the  CCSD  model  was  to  in  a  double 
zeta  basis.  A  potential  curve  from  0.9  a.u.  through  10  a.u.  was  calculated 
with  both  the  CCSO  model  and  SDCI  (i.e.  full  Cl),  which  of  course,  gave 
the  same  answer.  Furthermore,  the  CCSD  energy  was  invariant  to  the  orbital 
transformations  -  including  UHF  orbitals-  which  we  applied  to  test  the  program. 

The  second  application,  also  made  to  test  the  correctness  of  the 
program,  was  2H2  at  a  100  a.u.  intermolecular  separation.  In  this  problem, 
the  CCSD  energy  was  exactly  twice  the  SDCI  energy  for  a  single  H2  as 
size  extensivity  requires.  Again,  the  correct  relationship  held  for  the 
range  of  H2  intramolecular  separations  while  keeping  the  intermolecular 
separation  of  the  H^'s  at  100  a.u. 

More  interesting  results  are  the  CCSD  energies  for  l-^O  and  the 

comparison  with  the  full  Cl  results  of  Saxe,  Schaefer  and  Handy  (Table  5lf.24^0ur 

CCSD  energy,  in  the  double  zeta  basis  at  their  geometry,  is  -76.156077  hartrees 

which  represents  an  energy  lowering  of  -.006062  hartrees  relative  to  their 

SD- Cl  energy  of  -76.150015  hartrees.  Our  CCSD  energy  is  .001789  hartrees 

above  their  full  Cl  energy  of  -76.157866  hartrees.  The  largest  part  of  the 

difference  between  CCSD  and  full  Cl  comes  from  triple  excitations,  which 

(25) 

have  been  calculated  to  contribute  -0.001364  hartrees.  The  CCSD  energy  accounts 
for  98.79%  of  the  total  correlation  energy  compared  to  94.67%  for  SC-CI. 

Thus,  the  CCSD  model  accounts  for  two-thirds  of  the  error  inherent  in  SD-CI. 

As  can  be  seen  from  Table  5,  the  CCSD  +  T(4)  model  energy  is  .00017  hartrees 
(0.1  kcal/mole)  above  the  SDTQ-CI  energy.  In  this  case,  the  effect  of  Cl 
type  quadruple  excitations  are  handled  accurately  by  the  exponential  ansatz. 

* 


Equally  interesting  is  a  comparison  of  CCSD  with  CCSD-2,  a 


simplification  of  the  full  CCSD  in  which  only  the  linear  single  excitation 
2  (131 

terms  and  the  1/2T^  are  retained:  '  The  difference  between  CCSD  and  CCSD-2 
of  .00031  hartrees  is  the  energy  raising  caused  by  disconnected  single 
contributions  to  Cl  type  triple  and  quadruple  excitations. 

One  of  the  advantages  of  an  infinite-order  model  like  CCSD  over  a 
finite-order  perturbation  approximation  occurs  with  more  difficult  cases 
that  involve  some  quasidegeneracy!20^  Unlike  the  H ^0  example,  where  the 
Hartree-Fock  reference  determinant  has  a  coefficient  of  0.95  in  the  full 
Cl  expansion  with  all  other  coefficients  very  small,  a  quasidegenerate 
problem  usually  has  two  or  more  configurations  with  comparatively  large 
coefficients,  which  might  recommend  a  mul Preference  approach.  However, 
the  infinite-order  CCD  model  relative  to  only  a  single  reference  function 

(12) 

has  been  shown  to  often  describe  even  highly  degenerate  problems  reliably. 
Consequently,  to  assess  the  stability  of  the  CCSD  model  for  auasidegenerate  cases, 
a  problem  involving  the  insertion  of  Be  into  an  molecule  has  been  considered. 

Be  is  well-known  for  the  ouasidegeneracy  between  the  2s  and  2p 
2  2 

orbital  that  causes  the  Is  2p  configuration  to  be  important  in  the  Cl 

2  2 

expansion.  The  degree  of  importance  of  the  Is  2p  configuration  is  very 
much  a  function  of  the  choice  of  molecular  orbitals  with  MCSCF  orbitals 
weighting  it  heavily,  but  even  with  SCF  orbitals  from  a  large  basis  set,  the 

coefficient  is  still  about  O.l^^Also,  as  the  bond  is  broken,  the 
2  2 

lag  and  loy  configurations  become  equally  important.  In  addition,  to 
these  elements,  as  Be  is  inserted  perpendicularly  (Cjy)  into  H^,  there  is 
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a  promotion  from  Be(2s2)  to  Be(2p2)  near  the  critical  geometry,  causing 

2  2  2 

the  principal  configuration  for  BeHg  to  change  from  1  a -j  2a^  3a-|  to 
la2  2a2  lb2.  Since  in  a  single  reference  model,  one  of  these  configurations 
must  be  treated  in  the  complementary  space,  such  quasidegeneracy  effects 
should  severely  tax  the  ability  of  a  single  reference  model  to  describe 
this  insertion  reaction.  Also,  both  RHF  functions  are  unstable  since  a 
lower  UHF  solution  exists.  Just  as  in  full  Cl,  however,  even  a  poor  choice 
of  single  reference  function  might  be  used  to  generate  the  configurations, 
but  the  proper  weight  factors  would  be  obtained  via  the  diagonal ization  steD 
CCSD  potentially  has  the  same  flexibility,  although  one  must  distinguish 
between  a  method  containing  the  correct  solution  and  the  practical  problem 
of  extracting  that  solution  numerically.  In  all  applications  reported 
here  the  reduced  linear  equation  method  and  Pade  apDroximants^  ^  are  used 
to  obtain  the  solutions. 

The  basis  set  for  BeHg  is  given  in  Table  6,  and  we  present  the 
SCF,  full  Cl,  CCSD  and  various  lower  order  MBPT  results  in  Table  7. 

The  geometry, Be(0. ,  0. ,  3.0  a.u. )  and  H(0.  +1.16  a.u. ,  0. ) ,  is  near  the 

2  2  2 

point  of  crossover  when  the  principal  configuration  la^  2a-j  3a ^  would  be 
2  2  2 

superseded  by  la^  2a ^  lb^.  The  full  Cl  coefficients  using  SCF  orbitals 
2  2  2 

obtained  from  la-j  2a^  lb^  at  this  geometry  are  0.823  and  -.294  ,  respectively. 
When  Be  is  moved  to  2.75  a.u.  from  F^the  respective  coefficients  change 
to  -.560  and  0.724. 

Despite  the  large  amount  of  potential  degeneracy  in  this  system 
the  agreement  between  the  single  reference  CCSD  model  and  the  full  Cl  is 
exceptional.  Unlike  the  example  of  1^0  where  fourth-order  perturbation 
theory  (SDT0-MBPT(4)  is  only  0.6  kcal/mole  higher  than  the  full  Cl,  and 


the  other  fourth-order  approximations  are  in  the  vicinity  of  2  kcal/mole 
of  the  full  Cl,  the  error  is  6-7  kcal/mole  for  fourth-order  approximations 
for  BeH^. 

An  even  more  extreme  text  example  is  to  choose  to  use  the  less  important 
2  2  2 

la-j  2a-|  lb^  configuration  as  the  reference  determinant.  Its  coefficient 
of  0.294  corresponds  to  a  weight  of  less  than  10%  of  the  full  Cl  wavefunction. 
These  results  are  listed  in  Table  8.  Notice  that  the  SCF  result  is  much 
higher  and,  in  fact,  seems  to  be  approximating  the  second  root  of  the 
full  Cl.  The  perturbation  approximations  tend  to  cluster  in  the  vicinity 
of  the  second  eigenvalue  although  they  do  go  below  the  correct  answer, 
assuming  some  kind  of  average  value  between  the  two  eigenvalues.  However, 
the  CCSD  model  appears  to  overcome  the  comparatively  poor  starting  point 
to  a  large  degree,  getting  within  5.3  kcal/mole  of  the  full  Cl.  This  occurs 
despite  the  enormous  weight  of  the  la^  2a-j  3a^  configuration  (2.16  intermediately 
normalized)  in  the  CCSD  wavefunction.  BeH^  has  an  approximately  separated 
pair  structure  and  the  basis  set  is  small,  but  this  example  still  illustrates 
the  large  degree  of  flexibility  inherent  in  the  CCSD  model.  Combining  CCSD's 
stability,  with  its  efficiency,  size-extensivity,  and  its  equivalence  to 
the  full  Cl  for  the  chemically  pertinent  problem  of  a  group  of  separated 
electron  pairs,  appears  to  make  CCSD  a  very  attractive  model  for 
numerous  applications. 


n 


i 


APPENDIX  A 

The  Cl  Eigenvalue  Problem  in  Energy  Independent  Form 
Consider  the  problem  of  finding  an  eigenvalue  satisfying 

<H  -  EPS  ■  S 


(A-l) 


where  H  is  a  symmetric  matrix  and  is  the  eigenvector  corresponding  to  the 
energy  eigenvalue  E.  One  technique  for  solving  secular  equations  is  the 
partitioning  method.  Here  we  partition  equation  (A-l)  and  renormalize  it 
so  that 
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(A-2) 


which  is  equivalent  to  two  equations 


Equation  (A-6)  is  independent  of  E  and  equivalent  to  equation  (A-2).  ) 

The  energy  corresponding  to  the  t  which  satisfies  (A-6)  is  given  by  (A-5).  ; 

) 

The  quadratic  dependence  manifest  in  (A-6)  is  different  from  that  present  I 

in  the  coupled  cluster  equation,  since  in  (A-6)  there  is  a  scalar  product, 
while  in  the  coupled  cluster  equations  a  true  tensor  product,  t  x  t, 
appears. 

Equation  (A-6)  also  has  an  additional  interesting  feature  in  that 


the  linearized  coupled  cluster  method  can  be  derived  by  deleting  -(a  t_]_)t 
from  (A-6). 
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TABLE  2.  CCSD  equation  satisfied  by  the  double 
excitation  coefficient  t^jj 
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TABLE  3.  DEFINITIONS  OF  QUANTITIES  IN  TABLE  2 
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TABLE  4.  A  FORTRAN  SUBROUTINE  FOR  PERFORMING 
A  SPARSE  SCALAR  PRODUCT 

FUNCTION  PAKPRD  (A,  LA,  B,  LB,  N) 


PAKPRD  PERFORMS  A  PACKED  SCALAR  PRODUCT 
BETWEEN  VECTOR  A  AND  VECTOR  B 

A  :  AN  ARRAY  OF  NUMBERS  WITH  SKIP  INDICES 

LA  :  THE  LENGTH  OF  A 

B  :  AN  ARRAY  OF  NUMBERS  WITH  SKIP  INDICES 

LB  :  THE  LENGTH  OF  B 

N  :  THE  LENGTH  OF  A  AND  B  IF  THEY  WERE  UNPACKED 


PAKPRD  =  0. 

IA  =  0 

IB  s  0 

NA  =  0 

NB  *  0 

40  IA  ■  IA  +  1 

NA  -  NA  +  (A( IA)  -AND.  225) 

30  IB  =  IB  +  1 

NB  =  NB  +  ( B ( IB )  .AND.  255) 

50  IF  (NA-NB)  10,  20,  30 
20  PAKPRD  =  PAKPRD  +  A  (IA)*B(IB) 
IF  (NA  .LT,  N)  GO  TO  40 
IF  (NB  .GT.  N)  CALL  BOMB  (NB) 
IF  (IB  .GT,  LB)  CALL  BOMB  (IB) 
IF  (IA  .GT.  LA)  CALL  BOMB  (IA) 
RETURN 

10  IA  =  IA  +  1 

NA  =  NA  +  (A( IA)  .AND.  255) 

GO  TO  50 


TABLE  5.  COMPARISON  OF  MBPT/CCM  RESULTS  WITH  FULL  Cl 
H20  IN  14  CGTO  BASIS  SET  FROM  SAXE,  SCHAEFER 
AND  HANDY3 


Model 

Configurations 

Correlation 

Energy 

AE(FCI) 

Kcal/mole 

SD-CI 

361 

-0.14018 

4.9 

SDTQ-CI 

17,678 

-0.14777 

0.2 

FULL  Cl 

256,743 

-0.14803 

0.0 

D-MBPT(2) 

-0.13948 

5.4 

D-MBPT(3) 

-0.14087 

4.5 

D-MBPT(4) 

-0.14392 

2.6 

DQ-MBPT(4) 

-0.14476 

2.1 

SDQ-MBPT(4) 

-0.14565 

1.5 

SDTQ-MBPT(4)b 

-0.14704 

0.6 

CCD 

-0.14544 

1.6 

CCD-Orbital 

Optimized 

-0.14622 

1.3 

CCD  +  ST(4) 

-0.14771 

0.2 

CCSD-1  or  2 

-0.14655 

0.9 

CCSD 

-0.14624 

1.2 

CCSD  +  T(4) 

-0.14760 

0.3 

3  Reference  (24). 


b  Reference  (25). 


CONTRACTED  GAUSSIAN  BASIS  USED  FOR  THE  TEN 
ORBITAL  BeH2  MODEL  PROBLEM 


Contraction 


Exponent 

Coefficient 

Is 

1267.07 

.001940 

190.356 

.014786 

43.2959 

.071795 

12.1442 

.236348 

3.80923 

.471763 

1.26847 

.355183 

ls‘ 

5.693880 

-0.028876 

1.555630 

-0.177565 

0.171855 

1.071630 

Is" 

0.057181 

1.000 

2p 

5.693880 

.004836 

1.555630 

.144045 

0.171855 

.949692 

Is 

19.2406 

.032828 

2.8992 

.231208 

0.6534 

.817238 

Is' 

0.17760 

1.000 

TABLE  7.  COMPARISON  OF  BeH2  ENERGIES  CALCULATED  NEAR  THE 
TRANSITION  STATE  GEOMETRY  BE(0.,  0.,  3.  a.u.) 
H(0.,+1. 16a.  u.,0.)  REFERENCE  CONFIGURATION  a*a^ 


Model 

Configurations 

Total 

Energy 

aE(FCI) 

kcal/mole 

SCF 

1 

-15.53647 

55.5 

FCI 

1574 

-15.62496 

0.0 

D-MBPT(2) 

- 

-15.58485 

25.2 

D-MBPT(3) 

-15.60460 

12.8 

D-MBPT(4) 

-15.61437 

6.6 

SD-MBPT(4) 

-15.61485 

6.3 

DQ-MBPT(4) 

-15.61331 

7.3 

SDQ-MBPT(4) 

-15.61378 

7.0 

CCSD-2 

-15.62709 

1.3 

CCSD 

-15.62418 

0.5 

*  CCSD  expansion  coefficient  for  a^  b 2  =  -0.24 
using  intermediate  normalization. 


TABLE  8. 


COMPARISON  OF  BeH2  ENERGIES  CALCULATED  AT 
Be  (0. ,  0.  ,3a  .u. ) H  (0.,  +  1.16a.u.,  0.)  WITH 
REFERENCE  CONFIGURATION  a^i>2* 


Model 

Configurations 

Total 

Energy 

aE(FCI) 

kcal/mole 

SCF 

1 

-15.47728 

92.7 

FCI  Root  1 

1574 

-15.62496 

0.0 

FCI  Root  2 

1574 

-15.53575 

56.0 

D-MBPT(2) 

-15.51987 

65.9 

D-MBPT(3) 

-15.53564 

56.0 

D-MBPT(4) 

-15.54422 

50.7 

SD-MBPT(4) 

-15.54495 

50.2 

DQ-MBPT(4) 

-15.54331 

51.2 

SDQ-MBPT(4) 

-15.54404 

50.8 

CCSD 

-15.61645 

5.3 

*  CCSD  expansion 

coefficient  for 

a^a^  =  2.16 

using  intermediate  normalization. 
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